libs/numeric/odeint/performance/generic_odeint_rk4_lorenz.cpp

   1 /*
   2  * Copyright 2011 Mario Mulansky
   3  * Copyright 2012 Karsten Ahnert
   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 <boost/array.hpp>
  12
  13 //#include <boost/numeric/odeint/stepper/explicit_generic_rk.hpp>
  14 #include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
  15 #include <boost/numeric/odeint/algebra/array_algebra.hpp>
  16
  17 #include "rk_performance_test_case.hpp"
  18
  19 #include "lorenz.hpp"
  20
  21 using namespace boost::numeric::odeint;
  22
  23 typedef boost::array< double , 3 > state_type;
  24
  25 /*
  26 typedef explicit_generic_rk< 4 , 4 , state_type , double , state_type , double , array_algebra > rk4_type;
  27
  28 typedef rk4_type::coef_a_type coef_a_type;
  29 typedef rk4_type::coef_b_type coef_b_type;
  30 typedef rk4_type::coef_c_type coef_c_type;
  31
  32 const boost::array< double , 1 > a1 = {{ 0.5 }};
  33 const boost::array< double , 2 > a2 = {{ 0.0 , 0.5 }};
  34 const boost::array< double , 3 > a3 = {{ 0.0 , 0.0 , 1.0 }};
  35
  36 const coef_a_type a = fusion::make_vector( a1 , a2 , a3 );
  37 const coef_b_type b = {{ 1.0/6 , 1.0/3 , 1.0/3 , 1.0/6 }};
  38 const coef_c_type c = {{ 0.0 , 0.5 , 0.5 , 1.0 }};
  39 */
  40
  41 typedef runge_kutta4< state_type , double , state_type , double , array_algebra > rk4_type;
  42
  43
  44 class rk4_wrapper
  45 {
  46
  47 public:
  48
  49     rk4_wrapper()
  50     //        : m_stepper( a , b , c )
  51     {}
  52
  53     void reset_init_cond()
  54     {
  55         m_x[0] = 10.0 * rand() / RAND_MAX;
  56         m_x[1] = 10.0 * rand() / RAND_MAX;
  57         m_x[2] = 10.0 * rand() / RAND_MAX;
  58         m_t = 0.0;
  59     }
  60
  61     inline void do_step( const double dt )
  62     {
  63         m_stepper.do_step( lorenz(), m_x , m_t , dt );
  64     }
  65
  66     double state( const size_t i ) const
  67     { return m_x[i]; }
  68
  69 private:
  70     state_type m_x;
  71     double m_t;
  72     rk4_type m_stepper;
  73 };
  74
  75
  76
  77 int main()
  78 {
  79     srand( 12312354 );
  80
  81     rk4_wrapper stepper;
  82
  83     run( stepper );
  84 }