ceph/src/boost/libs/numeric/odeint/examples/bulirsch_stoer.cpp

   1 /*
   2  * bulirsch_stoer.cpp
   3  *
   4  * Copyright 2011-2013 Mario Mulansky
   5  * Copyright 2011-2012 Karsten Ahnert
   6  *
   7  * Distributed under the Boost Software License, Version 1.0.
   8  * (See accompanying file LICENSE_1_0.txt or
   9  * copy at http://www.boost.org/LICENSE_1_0.txt)
  10  */
  11
  12 #include <iostream>
  13 #include <fstream>
  14 #define _USE_MATH_DEFINES
  15 #include <cmath>
  16
  17 #include <boost/array.hpp>
  18 #include <boost/ref.hpp>
  19
  20 #include <boost/numeric/odeint/config.hpp>
  21
  22 #include <boost/numeric/odeint.hpp>
  23 #include <boost/numeric/odeint/stepper/bulirsch_stoer.hpp>
  24 #include <boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp>
  25
  26 using namespace std;
  27 using namespace boost::numeric::odeint;
  28
  29 typedef boost::array< double , 1 > state_type;
  30
  31 /*
  32  * x' = ( - x*sin t  + 2 tan x ) y
  33  * with x( pi/6 ) = 2/sqrt(3) the analytic solution is 1/cos t
  34  */
  35
  36 void rhs( const state_type &x , state_type &dxdt , const double t )
  37 {
  38     dxdt[0] = ( - x[0] * sin( t ) + 2.0 * tan( t ) ) * x[0];
  39 }
  40
  41 void rhs2( const state_type &x , state_type &dxdt , const double t )
  42 {
  43     dxdt[0] = sin(t);
  44 }
  45
  46
  47 ofstream out;
  48
  49 void write_out( const state_type &x , const double t )
  50 {
  51     out << t << '\t' << x[0] << endl;
  52 }
  53
  54 int main()
  55 {
  56     bulirsch_stoer_dense_out< state_type > stepper( 1E-8 , 0.0 , 0.0 , 0.0 );
  57     bulirsch_stoer< state_type > stepper2( 1E-8 , 0.0 , 0.0 , 0.0 );
  58
  59     state_type x = {{ 2.0 / sqrt(3.0) }};
  60
  61     double t = M_PI/6.0;
  62     //double t = 0.0;
  63     double dt = 0.01;
  64     double t_end = M_PI/2.0 - 0.1;
  65     //double t_end = 100.0;
  66
  67     out.open( "bs.dat" );
  68     out.precision(16);
  69     integrate_const( stepper , rhs , x , t , t_end , dt , write_out );
  70     out.close();
  71
  72     x[0] = 2.0 / sqrt(3.0);
  73
  74     out.open( "bs2.dat" );
  75     out.precision(16);
  76     integrate_adaptive( stepper , rhs , x , t , t_end , dt , write_out );
  77     out.close();
  78
  79     x[0] = 2.0 / sqrt(3.0);
  80
  81     out.open( "bs3.dat" );
  82     out.precision(16);
  83     integrate_adaptive( stepper2 , rhs , x , t , t_end , dt , write_out );
  84     out.close();
  85
  86
  87     typedef runge_kutta_dopri5< state_type > dopri5_type;
  88     typedef controlled_runge_kutta< dopri5_type > controlled_dopri5_type;
  89     typedef dense_output_runge_kutta< controlled_dopri5_type > dense_output_dopri5_type;
  90
  91     dense_output_dopri5_type dopri5 = make_dense_output( 1E-9 , 1E-9 , dopri5_type() );
  92
  93     x[0] = 2.0 / sqrt(3.0);
  94
  95     out.open( "bs4.dat" );
  96     out.precision(16);
  97     integrate_adaptive( dopri5 , rhs , x , t , t_end , dt , write_out );
  98     out.close();
  99
 100 }