implicit_euler.cpp 2.1 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586
  1. /*
  2. [auto_generated]
  3. libs/numeric/odeint/test/implicit_euler.cpp
  4. [begin_description]
  5. This file tests the implicit Euler stepper.
  6. [end_description]
  7. Copyright 2010-2011 Mario Mulansky
  8. Copyright 2010-2012 Karsten Ahnert
  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_implicit_euler
  19. #include <boost/test/unit_test.hpp>
  20. #include <utility>
  21. #include <iostream>
  22. #include <boost/numeric/odeint/stepper/implicit_euler.hpp>
  23. //#include <boost/numeric/odeint/util/ublas_resize.hpp>
  24. #include <boost/numeric/ublas/vector.hpp>
  25. #include <boost/numeric/ublas/matrix.hpp>
  26. using namespace boost::unit_test;
  27. using namespace boost::numeric::odeint;
  28. typedef double value_type;
  29. typedef boost::numeric::ublas::vector< value_type > state_type;
  30. typedef boost::numeric::ublas::matrix< value_type > matrix_type;
  31. /* use functors, because functions don't work with msvc 10, I guess this is a bug */
  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 ) const
  43. {
  44. jacobi( 0 , 0 ) = 1;
  45. jacobi( 0 , 1 ) = 2;
  46. jacobi( 1 , 0 ) = 0;
  47. jacobi( 1 , 1 ) = 1;
  48. }
  49. };
  50. BOOST_AUTO_TEST_SUITE( implicit_euler_test )
  51. BOOST_AUTO_TEST_CASE( test_euler )
  52. {
  53. implicit_euler< value_type > stepper;
  54. state_type x( 2 );
  55. x(0) = 0.0; x(1) = 1.0;
  56. value_type eps = 1E-12;
  57. /* make_pair doesn't work with function pointers on msvc 10 */
  58. stepper.do_step( std::make_pair( sys() , jacobi() ) , x , 0.0 , 0.1 );
  59. using std::abs;
  60. // compare with analytic solution of above system
  61. BOOST_CHECK_MESSAGE( abs( x(0) - 20.0/81.0 ) < eps , x(0) - 20.0/81.0 );
  62. BOOST_CHECK_MESSAGE( abs( x(1) - 10.0/9.0 ) < eps , x(0) - 10.0/9.0 );
  63. }
  64. BOOST_AUTO_TEST_SUITE_END()