libs/numeric/odeint/test/generic_error_stepper.cpp

   1 /*
   2  [auto_generated]
   3  libs/numeric/odeint/test/generic_error_stepper.cpp
   4
   5  [begin_description]
   6  This file tests the generic error stepper.
   7  [end_description]
   8
   9  Copyright 2011 Mario Mulansky
  10  Copyright 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 // disable checked iterator warning for msvc
  19 #include <boost/config.hpp>
  20 #ifdef BOOST_MSVC
  21     #pragma warning(disable:4996)
  22 #endif
  23
  24 #define BOOST_TEST_MODULE odeint_generic_error_stepper
  25
  26 #include <iostream>
  27 #include <utility>
  28
  29 #include <boost/test/unit_test.hpp>
  30
  31 #include <boost/numeric/odeint/stepper/explicit_error_generic_rk.hpp>
  32 #include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp>
  33 #include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54_classic.hpp>
  34
  35 #include <boost/array.hpp>
  36
  37 using namespace boost::unit_test;
  38 using namespace boost::numeric::odeint;
  39
  40 namespace fusion = boost::fusion;
  41
  42 typedef double value_type;
  43 typedef boost::array< value_type , 2 > state_type;
  44
  45 void sys( const state_type &x , state_type &dxdt , const value_type &t )
  46 {
  47     dxdt[ 0 ] = x[ 0 ] + 2 * x[ 1 ];
  48     dxdt[ 1 ] = x[ 1 ];
  49 }
  50
  51 typedef explicit_error_generic_rk< 6 , 5 , 4 , 5 , state_type> error_rk_generic_type;
  52
  53 const boost::array< double , 1 > a1 = {{ 0.2 }};
  54 const boost::array< double , 2 > a2 = {{ 3.0/40.0 , 9.0/40 }};
  55 const boost::array< double , 3 > a3 = {{ 0.3 , -0.9 , 1.2 }};
  56 const boost::array< double , 4 > a4 = {{ -11.0/54.0 , 2.5 , -70.0/27.0 , 35.0/27.0 }};
  57 const boost::array< double , 5 > a5 = {{ 1631.0/55296.0 , 175.0/512.0 , 575.0/13824.0 , 44275.0/110592.0 , 253.0/4096.0 }};
  58
  59 const error_rk_generic_type::coef_a_type a = fusion::make_vector( a1 , a2 , a3 , a4 , a5 );
  60 const error_rk_generic_type::coef_b_type b = {{ 37.0/378.0 , 0.0 , 250.0/621.0 , 125.0/594.0 , 0.0 , 512.0/1771.0 }};
  61 const error_rk_generic_type::coef_b_type b2 = {{ b[0]-2825.0/27648.0 , b[1]-0.0 , b[2]-18575.0/48384.0 , b[3]-13525.0/55296.0 , b[4]-277.0/14336.0 , b[5]-0.25 }};
  62 const error_rk_generic_type::coef_c_type c = {{ 0.0 , 0.2 , 0.3 , 0.6 , 1.0 , 7.0/8 }};
  63
  64 typedef runge_kutta_cash_karp54< state_type > error_rk54_ck_generic_type;
  65 typedef runge_kutta_cash_karp54_classic< state_type > rk54_ck_type;
  66
  67 BOOST_AUTO_TEST_SUITE( generic_error_stepper_test )
  68
  69 BOOST_AUTO_TEST_CASE( test_generic_error_stepper )
  70 {
  71     //simultaneously test copying
  72     error_rk_generic_type rk_generic_( a , b , b2 , c );
  73     error_rk_generic_type rk_generic = rk_generic_;
  74
  75     error_rk54_ck_generic_type rk54_ck_generic_;
  76     error_rk54_ck_generic_type rk54_ck_generic = rk54_ck_generic_;
  77
  78     //std::cout << stepper;
  79
  80     rk54_ck_type rk54_ck_;
  81     rk54_ck_type rk54_ck = rk54_ck_;
  82
  83     typedef error_rk_generic_type::state_type state_type;
  84     typedef error_rk_generic_type::value_type stepper_value_type;
  85     typedef error_rk_generic_type::deriv_type deriv_type;
  86     typedef error_rk_generic_type::time_type time_type;
  87
  88     state_type x = {{ 0.0 , 1.0 }};
  89     state_type y = x;
  90     state_type z = x;
  91     state_type x_err , y_err , z_err;
  92
  93     rk_generic.do_step( sys , x , 0.0 , 0.1 , x_err );
  94
  95     rk54_ck_generic.do_step( sys , y , 0.0 , 0.1 , y_err);
  96
  97     rk54_ck.do_step( sys , z , 0.0 , 0.1 , z_err );
  98
  99     BOOST_CHECK_NE( 0.0 , x[0] );
 100     BOOST_CHECK_NE( 1.0 , x[1] );
 101     BOOST_CHECK_NE( 0.0 , x_err[0] );
 102     BOOST_CHECK_NE( 0.0 , x_err[1] );
 103     // compare with analytic solution of above system
 104     BOOST_CHECK_EQUAL( x[0] , y[0] );
 105     BOOST_CHECK_EQUAL( x[1] , y[1] );
 106     BOOST_CHECK_EQUAL( x[0] , z[0] );
 107     BOOST_CHECK_EQUAL( x[1] , z[1] );
 108
 109 }
 110
 111 BOOST_AUTO_TEST_SUITE_END()