libs/numeric/odeint/test/integrate.cpp

   1 /*
   2  [auto_generated]
   3  libs/numeric/odeint/test/integrate.cpp
   4
   5  [begin_description]
   6  This file tests the integrate function and its variants.
   7  [end_description]
   8
   9  Copyright 2011-2012 Mario Mulansky
  10  Copyright 2011-2012 Karsten Ahnert
  11
  12  Distributed under the Boost Software License, Version 1.0.
  13  (See accompanying file LICENSE_1_0.txt or
  14  copy at http://www.boost.org/LICENSE_1_0.txt)
  15  */
  16
  17
  18 #define BOOST_TEST_MODULE odeint_integrate_functions
  19
  20 #include <vector>
  21 #include <cmath>
  22 #include <iostream>
  23
  24 #include <boost/numeric/odeint/config.hpp>
  25
  26 #include <boost/array.hpp>
  27 #include <boost/ref.hpp>
  28 #include <boost/iterator/counting_iterator.hpp>
  29
  30 #include <boost/test/unit_test.hpp>
  31
  32 #include <boost/mpl/vector.hpp>
  33
  34 // nearly everything from odeint is used in these tests
  35 #ifndef ODEINT_INTEGRATE_ITERATOR
  36 #include <boost/numeric/odeint/integrate/integrate_const.hpp>
  37 #include <boost/numeric/odeint/integrate/integrate_adaptive.hpp>
  38 #include <boost/numeric/odeint/integrate/integrate_times.hpp>
  39 #include <boost/numeric/odeint/integrate/integrate_n_steps.hpp>
  40 #else
  41 #include <boost/numeric/odeint/iterator/integrate/integrate_const.hpp>
  42 #include <boost/numeric/odeint/iterator/integrate/integrate_adaptive.hpp>
  43 #include <boost/numeric/odeint/iterator/integrate/integrate_times.hpp>
  44 #include <boost/numeric/odeint/iterator/integrate/integrate_n_steps.hpp>
  45 #endif
  46 #include <boost/numeric/odeint/stepper/euler.hpp>
  47 #include <boost/numeric/odeint/stepper/modified_midpoint.hpp>
  48 #include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
  49 #include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp>
  50 #include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
  51 #include <boost/numeric/odeint/stepper/runge_kutta_fehlberg78.hpp>
  52 #include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
  53 #include <boost/numeric/odeint/stepper/bulirsch_stoer.hpp>
  54 #include <boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp>
  55 #include <boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp>
  56
  57 #include <boost/numeric/odeint/util/detail/less_with_sign.hpp>
  58
  59 using namespace boost::unit_test;
  60 using namespace boost::numeric::odeint;
  61 namespace mpl = boost::mpl;
  62
  63
  64 typedef double value_type;
  65 typedef std::vector< value_type > state_type;
  66
  67 void lorenz( const state_type &x , state_type &dxdt , const value_type t )
  68 {
  69     //const value_type sigma( 10.0 );
  70     const value_type R( 28.0 );
  71     const value_type b( value_type( 8.0 ) / value_type( 3.0 ) );
  72
  73     // first component trivial
  74     dxdt[0] = 1.0; //sigma * ( x[1] - x[0] );
  75     dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
  76     dxdt[2] = -b * x[2] + x[0] * x[1];
  77 }
  78
  79 struct push_back_time
  80 {
  81     std::vector< double >& m_times;
  82
  83     state_type& m_x;
  84
  85     push_back_time( std::vector< double > &times , state_type &x )
  86     :  m_times( times ) , m_x( x ) { }
  87
  88     void operator()( const state_type &x , double t )
  89     {
  90         m_times.push_back( t );
  91         boost::numeric::odeint::copy( x , m_x );
  92     }
  93 };
  94
  95 template< class Stepper >
  96 struct perform_integrate_const_test
  97 {
  98     void operator()( const value_type t_end , const value_type dt )
  99     {
 100         std::cout << "Testing integrate_const with " << typeid( Stepper ).name() << std::endl;
 101
 102         state_type x( 3 , 10.0 ) , x_end( 3 );
 103
 104         std::vector< value_type > times;
 105
 106         size_t steps = integrate_const( Stepper() , lorenz , x , 0.0 , t_end ,
 107                          dt , push_back_time( times , x_end ) );
 108
 109         std::cout << t_end << " (" << dt << "), " << steps << " , " << times.size() << " , " << 10.0+dt*steps << "=" << x_end[0] << std::endl;
 110
 111         BOOST_CHECK_EQUAL( static_cast<int>(times.size()) , static_cast<int>(floor(t_end/dt))+1 );
 112
 113         for( size_t i=0 ; i<times.size() ; ++i )
 114         {
 115             //std::cout << i << " , " << times[i] << " , " << static_cast< value_type >(i)*dt << std::endl;
 116             // check if observer was called at times 0,1,2,...
 117             BOOST_CHECK_SMALL( times[i] - static_cast< value_type >(i)*dt , (i+1) * 2E-16 );
 118         }
 119
 120         // check first, trivial, component
 121         BOOST_CHECK_SMALL( (10.0 + times[times.size()-1]) - x_end[0] , 1E-6 ); // precision of steppers: 1E-6
 122         //BOOST_CHECK_EQUAL( x[1] , x_end[1] );
 123         //BOOST_CHECK_EQUAL( x[2] , x_end[2] );
 124     }
 125 };
 126
 127 template< class Stepper >
 128 struct perform_integrate_adaptive_test
 129 {
 130     void operator()( const value_type t_end = 10.0 , const value_type dt = 0.03 )
 131     {
 132         std::cout << "Testing integrate_adaptive with " << typeid( Stepper ).name() << std::endl;
 133
 134         state_type x( 3 , 10.0 ) , x_end( 3 );
 135
 136         std::vector< value_type > times;
 137
 138         size_t steps = integrate_adaptive( Stepper() , lorenz , x , 0.0 , t_end ,
 139                                         dt , push_back_time( times , x_end ) );
 140
 141         std::cout << t_end << " , " << steps << " , " << times.size() << " , " << dt << " , " << 10.0+t_end << "=" << x_end[0] << std::endl;
 142
 143         BOOST_CHECK_EQUAL( times.size() , steps+1 );
 144
 145         BOOST_CHECK_SMALL( times[0] - 0.0 , 2E-16 );
 146         BOOST_CHECK_SMALL( times[times.size()-1] - t_end , times.size() * 2E-16 );
 147
 148         // check first, trivial, component
 149         BOOST_CHECK_SMALL( (10.0 + t_end) - x_end[0] , 1E-6 ); // precision of steppers: 1E-6
 150 //        BOOST_CHECK_EQUAL( x[1] , x_end[1] );
 151 //        BOOST_CHECK_EQUAL( x[2] , x_end[2] );
 152     }
 153 };
 154
 155
 156 template< class Stepper >
 157 struct perform_integrate_times_test
 158 {
 159     void operator()( const int n = 10 , const int dn=1 , const value_type dt = 0.03 )
 160     {
 161         std::cout << "Testing integrate_times with " << typeid( Stepper ).name() << std::endl;
 162
 163         state_type x( 3 ) , x_end( 3 );
 164         x[0] = x[1] = x[2] = 10.0;
 165
 166         std::vector< double > times;
 167
 168         std::vector< double > obs_times( abs(n) );
 169         for( int i=0 ; boost::numeric::odeint::detail::less_with_sign( static_cast<double>(i) ,
 170                        static_cast<double>(obs_times.size()) ,
 171                        dt ) ; i+=dn )
 172         {
 173             obs_times[i] = i;
 174         }
 175         // simple stepper
 176         integrate_times( Stepper() , lorenz , x , obs_times.begin() , obs_times.end() ,
 177                     dt , push_back_time( times , x_end ) );
 178
 179         BOOST_CHECK_EQUAL( static_cast<int>(times.size()) , abs(n) );
 180
 181         for( size_t i=0 ; i<times.size() ; ++i )
 182             // check if observer was called at times 0,1,2,...
 183             BOOST_CHECK_EQUAL( times[i] , static_cast<double>(i) );
 184
 185         // check first, trivial, component
 186         BOOST_CHECK_SMALL( (10.0 + 1.0*times[times.size()-1]) - x_end[0] , 1E-6 ); // precision of steppers: 1E-6
 187 //        BOOST_CHECK_EQUAL( x[1] , x_end[1] );
 188 //        BOOST_CHECK_EQUAL( x[2] , x_end[2] );
 189     }
 190 };
 191
 192 template< class Stepper >
 193 struct perform_integrate_n_steps_test
 194 {
 195     void operator()( const int n = 200 , const value_type dt = 0.01 )
 196     {
 197         std::cout << "Testing integrate_n_steps with " << typeid( Stepper ).name() << std::endl;
 198
 199         state_type x( 3 ) , x_end( 3 );
 200         x[0] = x[1] = x[2] = 10.0;
 201
 202         std::vector< double > times;
 203
 204         // simple stepper
 205         value_type end_time = integrate_n_steps( Stepper() , lorenz , x , 0.0 , dt , n , push_back_time( times , x_end ) );
 206
 207         BOOST_CHECK_SMALL( end_time - n*dt , 2E-16 );
 208         BOOST_CHECK_EQUAL( static_cast<int>(times.size()) , n+1 );
 209
 210         for( size_t i=0 ; i<times.size() ; ++i )
 211             // check if observer was called at times 0,1,2,...
 212             BOOST_CHECK_SMALL( times[i] - static_cast< value_type >(i)*dt , 2E-16 );
 213
 214         // check first, trivial, component
 215         BOOST_CHECK_SMALL( (10.0 + end_time) - x_end[0] , 1E-6 ); // precision of steppers: 1E-6
 216 //        BOOST_CHECK_EQUAL( x[1] , x_end[1] );
 217 //        BOOST_CHECK_EQUAL( x[2] , x_end[2] );
 218
 219     }
 220 };
 221
 222
 223
 224 class stepper_methods : public mpl::vector<
 225     euler< state_type > ,
 226     modified_midpoint< state_type > ,
 227     runge_kutta4< state_type > ,
 228     runge_kutta_cash_karp54< state_type > ,
 229     runge_kutta_dopri5< state_type > ,
 230     runge_kutta_fehlberg78< state_type > ,
 231     controlled_runge_kutta< runge_kutta_cash_karp54< state_type > > ,
 232     controlled_runge_kutta< runge_kutta_dopri5< state_type > > ,
 233     controlled_runge_kutta< runge_kutta_fehlberg78< state_type > > ,
 234     bulirsch_stoer< state_type > ,
 235     dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > >
 236     //bulirsch_stoer_dense_out< state_type >
 237 > { };
 238
 239
 240
 241 BOOST_AUTO_TEST_SUITE( integrate_test )
 242
 243 BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_const_test_case , Stepper, stepper_methods )
 244 {
 245     perform_integrate_const_test< Stepper > tester;
 246     tester( 1.005 , 0.01 );
 247     tester( 1.0 , 0.01 );
 248     tester( 1.1 , 0.01 );
 249     tester( -1.005 , -0.01 );
 250 }
 251
 252
 253 BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_adaptive_test_case , Stepper, stepper_methods )
 254 {
 255     perform_integrate_adaptive_test< Stepper > tester;
 256     tester( 1.005 , 0.01 );
 257     tester( 1.0 , 0.01 );
 258     tester( 1.1 , 0.01 );
 259     tester( -1.005 , -0.01 );
 260 }
 261
 262
 263 BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_times_test_case , Stepper, stepper_methods )
 264 {
 265     perform_integrate_times_test< Stepper > tester;
 266     tester();
 267     //tester( -10 , -0.01 );
 268 }
 269
 270 BOOST_AUTO_TEST_CASE_TEMPLATE( integrate_n_steps_test_case , Stepper, stepper_methods )
 271 {
 272     perform_integrate_n_steps_test< Stepper > tester;
 273     tester();
 274     tester( 200 , 0.01 );
 275     tester( 200 , 0.01 );
 276     tester( 200 , 0.01 );
 277     tester( 200 , -0.01 );
 278 }
 279
 280 BOOST_AUTO_TEST_SUITE_END()