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
rhs(const state_type & x,state_type & dxdt,const double t)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
rhs2(const state_type & x,state_type & dxdt,const double t)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
write_out(const state_type & x,const double t)49 void write_out( const state_type &x , const double t )
50 {
51 out << t << '\t' << x[0] << endl;
52 }
53
main()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 }
101