rosenbrock4.cpp 5.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173
  1. /*
  2. [auto_generated]
  3. libs/numeric/odeint/test/rosenbrock4.cpp
  4. [begin_description]
  5. This file tests the Rosenbrock 4 stepper and its controller and dense output stepper.
  6. [end_description]
  7. Copyright 2011-2012 Karsten Ahnert
  8. Copyright 2011 Mario Mulansky
  9. Distributed under the Boost Software License, Version 1.0.
  10. (See accompanying file LICENSE_1_0.txt or
  11. copy at http://www.boost.org/LICENSE_1_0.txt)
  12. */
  13. // disable checked iterator warning for msvc
  14. #include <boost/config.hpp>
  15. #ifdef BOOST_MSVC
  16. #pragma warning(disable:4996)
  17. #endif
  18. #define BOOST_TEST_MODULE odeint_rosenbrock4
  19. #include <utility>
  20. #include <iostream>
  21. #include <boost/test/unit_test.hpp>
  22. #include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
  23. #include <boost/numeric/odeint/stepper/rosenbrock4_controller.hpp>
  24. #include <boost/numeric/odeint/stepper/rosenbrock4_dense_output.hpp>
  25. #include <boost/numeric/ublas/vector.hpp>
  26. #include <boost/numeric/ublas/matrix.hpp>
  27. using namespace boost::unit_test;
  28. using namespace boost::numeric::odeint;
  29. typedef double value_type;
  30. typedef boost::numeric::ublas::vector< value_type > state_type;
  31. typedef boost::numeric::ublas::matrix< value_type > matrix_type;
  32. struct sys
  33. {
  34. void operator()( const state_type &x , state_type &dxdt , const value_type &t ) const
  35. {
  36. dxdt( 0 ) = x( 0 ) + 2 * x( 1 );
  37. dxdt( 1 ) = x( 1 );
  38. }
  39. };
  40. struct jacobi
  41. {
  42. void operator()( const state_type &x , matrix_type &jacobi , const value_type &t , state_type &dfdt ) const
  43. {
  44. jacobi( 0 , 0 ) = 1;
  45. jacobi( 0 , 1 ) = 2;
  46. jacobi( 1 , 0 ) = 0;
  47. jacobi( 1 , 1 ) = 1;
  48. dfdt( 0 ) = 0.0;
  49. dfdt( 1 ) = 0.0;
  50. }
  51. };
  52. BOOST_AUTO_TEST_SUITE( rosenbrock4_test )
  53. BOOST_AUTO_TEST_CASE( test_rosenbrock4_stepper )
  54. {
  55. typedef rosenbrock4< value_type > stepper_type;
  56. stepper_type stepper;
  57. typedef stepper_type::state_type state_type;
  58. typedef stepper_type::value_type stepper_value_type;
  59. typedef stepper_type::deriv_type deriv_type;
  60. typedef stepper_type::time_type time_type;
  61. state_type x( 2 ) , xerr( 2 );
  62. x(0) = 0.0; x(1) = 1.0;
  63. stepper.do_step( std::make_pair( sys() , jacobi() ) , x , 0.0 , 0.1 , xerr );
  64. stepper.do_step( std::make_pair( sys() , jacobi() ) , x , 0.0 , 0.1 );
  65. // using std::abs;
  66. // value_type eps = 1E-12;
  67. //
  68. // // compare with analytic solution of above system
  69. // BOOST_CHECK_SMALL( abs( x(0) - 20.0/81.0 ) , eps );
  70. // BOOST_CHECK_SMALL( abs( x(1) - 10.0/9.0 ) , eps );
  71. }
  72. BOOST_AUTO_TEST_CASE( test_rosenbrock4_controller )
  73. {
  74. typedef rosenbrock4_controller< rosenbrock4< value_type > > stepper_type;
  75. stepper_type stepper;
  76. typedef stepper_type::state_type state_type;
  77. typedef stepper_type::value_type stepper_value_type;
  78. typedef stepper_type::deriv_type deriv_type;
  79. typedef stepper_type::time_type time_type;
  80. state_type x( 2 );
  81. x( 0 ) = 0.0 ; x(1) = 1.0;
  82. value_type t = 0.0 , dt = 0.01;
  83. stepper.try_step( std::make_pair( sys() , jacobi() ) , x , t , dt );
  84. }
  85. BOOST_AUTO_TEST_CASE( test_rosenbrock4_dense_output )
  86. {
  87. typedef rosenbrock4_dense_output< rosenbrock4_controller< rosenbrock4< value_type > > > stepper_type;
  88. typedef rosenbrock4_controller< rosenbrock4< value_type > > controlled_stepper_type;
  89. controlled_stepper_type c_stepper;
  90. stepper_type stepper( c_stepper );
  91. typedef stepper_type::state_type state_type;
  92. typedef stepper_type::value_type stepper_value_type;
  93. typedef stepper_type::deriv_type deriv_type;
  94. typedef stepper_type::time_type time_type;
  95. state_type x( 2 );
  96. x( 0 ) = 0.0 ; x(1) = 1.0;
  97. stepper.initialize( x , 0.0 , 0.1 );
  98. std::pair< value_type , value_type > tr = stepper.do_step( std::make_pair( sys() , jacobi() ) );
  99. stepper.calc_state( 0.5 * ( tr.first + tr.second ) , x );
  100. }
  101. class rosenbrock4_controller_max_dt_adaptable : public rosenbrock4_controller< rosenbrock4< value_type > >
  102. {
  103. public:
  104. void set_max_dt(value_type max_dt)
  105. {
  106. m_max_dt = max_dt;
  107. }
  108. };
  109. BOOST_AUTO_TEST_CASE( test_rosenbrock4_dense_output_ref )
  110. {
  111. typedef rosenbrock4_dense_output< boost::reference_wrapper< rosenbrock4_controller_max_dt_adaptable > > stepper_type;
  112. rosenbrock4_controller_max_dt_adaptable c_stepper;
  113. stepper_type stepper( boost::ref( c_stepper ) );
  114. typedef stepper_type::state_type state_type;
  115. typedef stepper_type::value_type stepper_value_type;
  116. typedef stepper_type::deriv_type deriv_type;
  117. typedef stepper_type::time_type time_type;
  118. state_type x( 2 );
  119. x( 0 ) = 0.0 ; x(1) = 1.0;
  120. stepper.initialize( x , 0.0 , 0.1 );
  121. std::pair< value_type , value_type > tr = stepper.do_step( std::make_pair( sys() , jacobi() ) );
  122. stepper.calc_state( 0.5 * ( tr.first + tr.second ) , x );
  123. // adapt the maximal step size to a small value (smaller than the step size) and make a step
  124. const double max_dt = 1e-8;
  125. c_stepper.set_max_dt(max_dt);
  126. stepper.do_step( std::make_pair( sys() , jacobi() ) );
  127. // check if the step was done with the requested maximal step size
  128. BOOST_CHECK_CLOSE(max_dt, stepper.current_time_step(), 1e-14);
  129. }
  130. BOOST_AUTO_TEST_CASE( test_rosenbrock4_copy_dense_output )
  131. {
  132. typedef rosenbrock4_controller< rosenbrock4< value_type > > controlled_stepper_type;
  133. typedef rosenbrock4_dense_output< controlled_stepper_type > stepper_type;
  134. controlled_stepper_type c_stepper;
  135. stepper_type stepper( c_stepper );
  136. stepper_type stepper2( stepper );
  137. }
  138. BOOST_AUTO_TEST_SUITE_END()