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

   1 /* Boost numeric test of the adams-bashforth steppers test file
   2
   3  Copyright 2013 Karsten Ahnert
   4  Copyright 2013-2015 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_adaptive_adams_bashforth_moulton
  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/mpl/vector.hpp>
  27
  28 #include <boost/numeric/odeint.hpp>
  29
  30 using namespace boost::unit_test;
  31 using namespace boost::numeric::odeint;
  32 namespace mpl = boost::mpl;
  33
  34 typedef double value_type;
  35
  36 typedef boost::array< double , 2 > state_type;
  37 typedef runge_kutta_fehlberg78<state_type> initializing_stepper;
  38
  39 // harmonic oscillator, analytic solution x[0] = sin( t )
  40 struct osc
  41 {
  42     void operator()( const state_type &x , state_type &dxdt , const double t ) const
  43     {
  44         dxdt[0] = x[1];
  45         dxdt[1] = -x[0];
  46     }
  47 };
  48
  49 BOOST_AUTO_TEST_SUITE( numeric_adaptive_adams_bashforth_moulton_test )
  50
  51
  52 /* generic test for all adams bashforth steppers */
  53 template< class Stepper >
  54 struct perform_adaptive_adams_bashforth_moulton_test
  55 {
  56     void operator()( void )
  57     {
  58         Stepper stepper;
  59         initializing_stepper init_stepper;
  60
  61         const int o = stepper.order()+1; //order of the error is order of approximation + 1
  62
  63         const state_type x0 = {{ 0.0 , 1.0 }};
  64         state_type x1 = x0;
  65         double t = 0.0;
  66         double dt = 0.25;
  67         // initialization, does a number of steps to self-start the stepper with a small stepsize
  68         stepper.initialize( init_stepper, osc() , x1 , t ,  dt);
  69         double A = std::sqrt( x1[0]*x1[0] + x1[1]*x1[1] );
  70         double phi = std::asin(x1[0]/A) - t;
  71
  72         // now we do the actual step
  73         stepper.do_step( osc() , x1 , t , dt );
  74         // only examine the error of the adams-bashforth step, not the initialization
  75         const double f = 2.0 * std::abs( A*sin(t+dt+phi) - x1[0] ) / std::pow( dt , o ); // upper bound
  76
  77         std::cout << o << " , " << f << std::endl;
  78
  79         /* as long as we have errors above machine precision */
  80         while( f*std::pow( dt , o ) > 1E-16 )
  81         {
  82             x1 = x0;
  83             t = 0.0;
  84             stepper.initialize( init_stepper, osc() , x1 , t , dt );
  85             A = std::sqrt( x1[0]*x1[0] + x1[1]*x1[1] );
  86             phi = std::asin(x1[0]/A) - t;
  87             // now we do the actual step
  88             stepper.do_step( osc() , x1 , t , dt );
  89             stepper.reset();
  90             // only examine the error of the adams-bashforth step, not the initialization
  91             std::cout << "Testing dt=" << dt << " , " << std::abs( A*sin(t+dt+phi) - x1[0] ) << std::endl;
  92             BOOST_CHECK_LT( std::abs( A*sin(t+dt+phi) - x1[0] ) , f*std::pow( dt , o ) );
  93             dt *= 0.5;
  94         }
  95     }
  96 };
  97
  98 typedef mpl::vector<
  99     adaptive_adams_bashforth_moulton< 2 , state_type > ,
 100     adaptive_adams_bashforth_moulton< 3 , state_type > ,
 101     adaptive_adams_bashforth_moulton< 4 , state_type > ,
 102     adaptive_adams_bashforth_moulton< 5 , state_type > ,
 103     adaptive_adams_bashforth_moulton< 6 , state_type > ,
 104     adaptive_adams_bashforth_moulton< 7 , state_type >
 105     > adaptive_adams_bashforth_moulton_steppers;
 106
 107 BOOST_AUTO_TEST_CASE_TEMPLATE( adaptive_adams_bashforth_moulton_test , Stepper, adaptive_adams_bashforth_moulton_steppers )
 108 {
 109     perform_adaptive_adams_bashforth_moulton_test< Stepper > tester;
 110     tester();
 111 }
 112
 113 BOOST_AUTO_TEST_SUITE_END()