ceph/src/boost/libs/numeric/odeint/test/numeric/rosenbrock.cpp

   1 /* Boost numeric test of the rosenbrock4 stepper test file
   2
   3  Copyright 2012 Karsten Ahnert
   4  Copyright 2012 Mario Mulansky
   5
   6  Distributed under the Boost Software License, Version 1.0.
   7  (See accompanying file LICENSE_1_0.txt or
   8  copy at http://www.boost.org/LICENSE_1_0.txt)
   9 */
  10
  11 // disable checked iterator warning for msvc
  12 #include <boost/config.hpp>
  13 #ifdef BOOST_MSVC
  14     #pragma warning(disable:4996)
  15 #endif
  16
  17 #define BOOST_TEST_MODULE numeric_rosenbrock
  18
  19 #include <iostream>
  20 #include <cmath>
  21
  22 #include <boost/array.hpp>
  23
  24 #include <boost/test/unit_test.hpp>
  25
  26 #include <boost/numeric/odeint/stepper/rosenbrock4.hpp>
  27 #include <boost/numeric/ublas/vector.hpp>
  28 #include <boost/numeric/ublas/matrix.hpp>
  29
  30 using namespace boost::unit_test;
  31 using namespace boost::numeric::odeint;
  32
  33 typedef double value_type;
  34 typedef boost::numeric::ublas::vector< value_type > state_type;
  35 typedef boost::numeric::ublas::matrix< value_type > matrix_type;
  36
  37 // harmonic oscillator, analytic solution x[0] = sin( t )
  38 struct sys
  39 {
  40     void operator()( const state_type &x , state_type &dxdt , const value_type &t ) const
  41     {
  42         dxdt( 0 ) = x( 1 );
  43         dxdt( 1 ) = -x( 0 );
  44     }
  45 };
  46
  47 struct jacobi
  48 {
  49     void operator()( const state_type &x , matrix_type &jacobi , const value_type &t , state_type &dfdt ) const
  50     {
  51         jacobi( 0 , 0 ) = 0;
  52         jacobi( 0 , 1 ) = 1;
  53         jacobi( 1 , 0 ) = -1;
  54         jacobi( 1 , 1 ) = 0;
  55         dfdt( 0 ) = 0.0;
  56         dfdt( 1 ) = 0.0;
  57     }
  58 };
  59
  60
  61 BOOST_AUTO_TEST_SUITE( numeric_rosenbrock4 )
  62
  63 BOOST_AUTO_TEST_CASE( rosenbrock4_numeric_test )
  64 {
  65     typedef rosenbrock4< value_type > stepper_type;
  66     stepper_type stepper;
  67
  68     const int o = stepper.order()+1;
  69
  70     state_type x0( 2 ) , x1( 2 );
  71     x0(0) = 0.0; x0(1) = 1.0;
  72
  73     double dt = 0.5;
  74
  75     stepper.do_step( std::make_pair( sys() , jacobi() ) , x0 , 0.0 , x1 , dt );
  76     const double f = 2.0 * std::abs( std::sin(dt) - x1(0) ) / std::pow( dt , o );
  77
  78     std::cout << o << " , " << f << std::endl;
  79
  80     while( f*std::pow( dt , o ) > 1E-16 )
  81     {
  82         stepper.do_step( std::make_pair( sys() , jacobi() ) , x0 , 0.0 , x1 , dt );
  83         std::cout << "Testing dt=" << dt << std::endl;
  84         BOOST_CHECK_SMALL( std::abs( std::sin(dt) - x1(0) ) , f*std::pow( dt , o ) );
  85         dt *= 0.5;
  86     }
  87 }
  88
  89 BOOST_AUTO_TEST_SUITE_END()