ceph/src/boost/libs/numeric/odeint/examples/multiprecision/cmp_precision.cpp

   1 /* Boost libs/numeric/odeint/examples/multiprecision/cmp_precision.cpp
   2
   3  Copyright 2013 Karsten Ahnert
   4  Copyright 2013 Mario Mulansky
   5
   6  example comparing double to multiprecision using Boost.Multiprecision
   7
   8  Distributed under the Boost Software License, Version 1.0.
   9 (See accompanying file LICENSE_1_0.txt or
  10  copy at http://www.boost.org/LICENSE_1_0.txt)
  11  */
  12
  13
  14 #include <iostream>
  15 #include <boost/numeric/odeint.hpp>
  16 #include <boost/multiprecision/cpp_dec_float.hpp>
  17
  18 using namespace std;
  19 using namespace boost::numeric::odeint;
  20
  21 typedef boost::multiprecision::cpp_dec_float_50 mp_50;
  22
  23 /* we solve the simple ODE x' = 3/(2t^2) + x/(2t)
  24  * with initial condition x(1) = 0.
  25  * Analytic solution is x(t) = sqrt(t) - 1/t
  26  */
  27
  28 void rhs_m( const mp_50 x , mp_50 &dxdt , const mp_50 t )
  29 {   // version for multiprecision
  30     dxdt = mp_50(3)/(mp_50(2)*t*t) + x/(mp_50(2)*t);
  31 }
  32
  33 void rhs_d( const double x , double &dxdt , const double t )
  34 {   // version for double precision
  35     dxdt = 3.0/(2.0*t*t) + x/(2.0*t);
  36 }
  37
  38 // state_type = mp_50 = deriv_type = time_type = mp_50
  39 typedef runge_kutta4< mp_50 , mp_50 , mp_50 , mp_50 , vector_space_algebra , default_operations , never_resizer > stepper_type_m;
  40
  41 typedef runge_kutta4< double , double , double , double , vector_space_algebra , default_operations , never_resizer > stepper_type_d;
  42
  43 int main()
  44 {
  45
  46     stepper_type_m stepper_m;
  47     stepper_type_d stepper_d;
  48
  49     mp_50 dt_m( 0.5 );
  50     double dt_d( 0.5 );
  51
  52     cout << "dt" << '\t' << "mp" << '\t' << "double" << endl;
  53
  54     while( dt_m > 1E-20 )
  55     {
  56
  57         mp_50 x_m = 0; //initial value x(1) = 0
  58         stepper_m.do_step( rhs_m , x_m , mp_50( 1 ) , dt_m );
  59         double x_d = 0;
  60         stepper_d.do_step( rhs_d , x_d , 1.0 , dt_d );
  61
  62         cout << dt_m << '\t';
  63         cout << abs((x_m - (sqrt(1+dt_m)-mp_50(1)/(1+dt_m)))/x_m) << '\t' ;
  64         cout << abs((x_d - (sqrt(1+dt_d)-mp_50(1)/(1+dt_d)))/x_d) << endl ;
  65         dt_m /= 2;
  66         dt_d /= 2;
  67     }
  68 }