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

   1 /* Boost numeric test for orders of quadrature formulas test file
   2
   3  Copyright 2015 Gregor de Cillia
   4  Copyright 2015 Mario Mulansky <mario.mulansky@gmx.net>
   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
  14 #ifdef BOOST_MSVC
  15     #pragma warning(disable:4996)
  16 #endif
  17
  18 #define BOOST_TEST_MODULE order_quadrature_formula
  19
  20 #include <iostream>
  21 #include <cmath>
  22
  23 #include "boost/format.hpp"
  24
  25 #include <boost/test/unit_test.hpp>
  26
  27 #include <boost/mpl/vector.hpp>
  28
  29 #include <boost/numeric/odeint.hpp>
  30
  31 #include <boost/numeric/ublas/vector.hpp>
  32
  33 #include <boost/format.hpp>
  34
  35 using namespace boost::unit_test;
  36 using namespace boost::numeric::odeint;
  37 namespace mpl = boost::mpl;
  38
  39 typedef double value_type;
  40 typedef value_type time_type;
  41 typedef value_type state_type;
  42
  43 BOOST_AUTO_TEST_SUITE( order_of_convergence_test )
  44
  45 /* defines the simple monomial f(t) = (p+1) * t^p.*/
  46 struct monomial
  47 {
  48     int power;
  49
  50     monomial(int p = 0) : power( p ){};
  51
  52     void operator()( const state_type &x , state_type &dxdt , const time_type t )
  53     {
  54         dxdt = ( 1.0 + power ) * pow( t, power );
  55     }
  56 };
  57
  58
  59 /* generic test for all steppers that support integrate_const */
  60 template< class Stepper >
  61 struct stepper_order_test
  62 {
  63     void operator()( int steps = 1 )
  64     {
  65         const int estimated_order = estimate_order( steps );
  66         const int defined_order = Stepper::order_value;
  67
  68         std::cout << boost::format( "%-20i%-20i\n" )
  69             % estimated_order %  defined_order;
  70
  71         BOOST_REQUIRE_EQUAL( estimated_order, defined_order );
  72     }
  73
  74     /*
  75     the order of the stepper is estimated by trying to solve the ODE
  76     x'(t) = (p+1) * t^p
  77     until the errors are too big to be justified by finite precision.
  78     the first value p for which the problem is *not* solved within the
  79     finite precision tolerance is the estimate for the order of the scheme.
  80      */
  81     int estimate_order( int steps )
  82     {
  83         const double dt = 1.0/steps;
  84         const double tolerance = steps*1E-15;
  85         int p;
  86         for( p = 0; true; p++ )
  87         {
  88             // begin with x'(t) = t^0 = 1
  89             //         => x (t) = t
  90             // then use   x'(t) = 2*t^1
  91             //         => x (t) = t^2
  92             // ...
  93             state_type x = 0.0;
  94
  95             double t = integrate_n_steps( Stepper(), monomial( p ), x, 0.0, dt,
  96                                           steps );
  97             if( fabs( x - pow( t, ( 1.0 + p ) ) ) > tolerance )
  98                 break;
  99         }
 100         // the smallest power p for which the test failed is the estimated order,
 101         // as the solution for this power is x(t) = t^{p+1}
 102         return p;
 103     }
 104 };
 105
 106
 107 typedef mpl::vector<
 108     euler< state_type > ,
 109     modified_midpoint< state_type > ,
 110     runge_kutta4< state_type > ,
 111     runge_kutta4_classic< state_type > ,
 112     runge_kutta_cash_karp54_classic< state_type > ,
 113     runge_kutta_cash_karp54< state_type > ,
 114     runge_kutta_dopri5< state_type > ,
 115     runge_kutta_fehlberg78< state_type >
 116     > runge_kutta_steppers;
 117
 118 typedef mpl::vector<
 119     adams_bashforth< 2, state_type, double, state_type, double,
 120                      vector_space_algebra, default_operations,
 121                      initially_resizer, runge_kutta_fehlberg78< state_type > >,
 122     adams_bashforth< 3, state_type, double, state_type, double,
 123                      vector_space_algebra, default_operations,
 124                      initially_resizer, runge_kutta_fehlberg78< state_type > >,
 125     adams_bashforth< 4, state_type, double, state_type, double,
 126                      vector_space_algebra, default_operations,
 127                      initially_resizer, runge_kutta_fehlberg78< state_type > >,
 128     adams_bashforth< 5, state_type, double, state_type, double,
 129                      vector_space_algebra, default_operations,
 130                      initially_resizer, runge_kutta_fehlberg78< state_type > >,
 131     adams_bashforth< 6, state_type, double, state_type, double,
 132                      vector_space_algebra, default_operations,
 133                      initially_resizer, runge_kutta_fehlberg78< state_type > >,
 134     adams_bashforth< 7, state_type, double, state_type, double,
 135                      vector_space_algebra, default_operations,
 136                      initially_resizer, runge_kutta_fehlberg78< state_type > >,
 137     adams_bashforth< 8, state_type, double, state_type, double,
 138                      vector_space_algebra, default_operations,
 139                      initially_resizer, runge_kutta_fehlberg78< state_type > >
 140     > ab_steppers;
 141
 142
 143 typedef mpl::vector<
 144     adams_bashforth_moulton< 2, state_type, double, state_type, double,
 145                              vector_space_algebra, default_operations,
 146                              initially_resizer,
 147                              runge_kutta_fehlberg78< state_type > >,
 148     adams_bashforth_moulton< 3, state_type, double, state_type, double,
 149                              vector_space_algebra, default_operations,
 150                              initially_resizer,
 151                              runge_kutta_fehlberg78< state_type > >,
 152     adams_bashforth_moulton< 4, state_type, double, state_type, double,
 153                              vector_space_algebra, default_operations,
 154                              initially_resizer,
 155                              runge_kutta_fehlberg78< state_type > >,
 156     adams_bashforth_moulton< 5, state_type, double, state_type, double,
 157                              vector_space_algebra, default_operations,
 158                              initially_resizer,
 159                              runge_kutta_fehlberg78< state_type > >,
 160     adams_bashforth_moulton< 6, state_type, double, state_type, double,
 161                              vector_space_algebra, default_operations,
 162                              initially_resizer,
 163                              runge_kutta_fehlberg78< state_type > >,
 164     adams_bashforth_moulton< 7, state_type, double, state_type, double,
 165                              vector_space_algebra, default_operations,
 166                              initially_resizer,
 167                              runge_kutta_fehlberg78< state_type > >,
 168     adams_bashforth_moulton< 8, state_type, double, state_type, double,
 169                              vector_space_algebra, default_operations,
 170                              initially_resizer,
 171                              runge_kutta_fehlberg78< state_type > >
 172     > abm_steppers;
 173
 174
 175 BOOST_AUTO_TEST_CASE_TEMPLATE( runge_kutta_test , Stepper, runge_kutta_steppers )
 176 {
 177     stepper_order_test< Stepper > tester;
 178     tester(10);
 179 }
 180
 181
 182 BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_test , Stepper, ab_steppers )
 183 {
 184     stepper_order_test< Stepper > tester;
 185     tester(16);
 186 }
 187
 188
 189 BOOST_AUTO_TEST_CASE_TEMPLATE( adams_bashforth_moultion_test , Stepper, abm_steppers )
 190 {
 191     stepper_order_test< Stepper > tester;
 192     tester(16);
 193 }
 194
 195 BOOST_AUTO_TEST_SUITE_END()