| /* |
| * bulirsch_stoer.cpp |
| * |
| * Copyright 2011-2013 Mario Mulansky |
| * Copyright 2011-2012 Karsten Ahnert |
| * |
| * Distributed under the Boost Software License, Version 1.0. |
| * (See accompanying file LICENSE_1_0.txt or |
| * copy at http://www.boost.org/LICENSE_1_0.txt) |
| */ |
| |
| #include <iostream> |
| #include <fstream> |
| #define _USE_MATH_DEFINES |
| #include <cmath> |
| |
| #include <boost/array.hpp> |
| #include <boost/ref.hpp> |
| |
| #include <boost/numeric/odeint/config.hpp> |
| |
| #include <boost/numeric/odeint.hpp> |
| #include <boost/numeric/odeint/stepper/bulirsch_stoer.hpp> |
| #include <boost/numeric/odeint/stepper/bulirsch_stoer_dense_out.hpp> |
| |
| using namespace std; |
| using namespace boost::numeric::odeint; |
| |
| typedef boost::array< double , 1 > state_type; |
| |
| /* |
| * x' = ( - x*sin t + 2 tan x ) y |
| * with x( pi/6 ) = 2/sqrt(3) the analytic solution is 1/cos t |
| */ |
| |
| void rhs( const state_type &x , state_type &dxdt , const double t ) |
| { |
| dxdt[0] = ( - x[0] * sin( t ) + 2.0 * tan( t ) ) * x[0]; |
| } |
| |
| void rhs2( const state_type &x , state_type &dxdt , const double t ) |
| { |
| dxdt[0] = sin(t); |
| } |
| |
| |
| ofstream out; |
| |
| void write_out( const state_type &x , const double t ) |
| { |
| out << t << '\t' << x[0] << endl; |
| } |
| |
| int main() |
| { |
| bulirsch_stoer_dense_out< state_type > stepper( 1E-8 , 0.0 , 0.0 , 0.0 ); |
| bulirsch_stoer< state_type > stepper2( 1E-8 , 0.0 , 0.0 , 0.0 ); |
| |
| state_type x = {{ 2.0 / sqrt(3.0) }}; |
| |
| double t = M_PI/6.0; |
| //double t = 0.0; |
| double dt = 0.01; |
| double t_end = M_PI/2.0 - 0.1; |
| //double t_end = 100.0; |
| |
| out.open( "bs.dat" ); |
| out.precision(16); |
| integrate_const( stepper , rhs , x , t , t_end , dt , write_out ); |
| out.close(); |
| |
| x[0] = 2.0 / sqrt(3.0); |
| |
| out.open( "bs2.dat" ); |
| out.precision(16); |
| integrate_adaptive( stepper , rhs , x , t , t_end , dt , write_out ); |
| out.close(); |
| |
| x[0] = 2.0 / sqrt(3.0); |
| |
| out.open( "bs3.dat" ); |
| out.precision(16); |
| integrate_adaptive( stepper2 , rhs , x , t , t_end , dt , write_out ); |
| out.close(); |
| |
| |
| typedef runge_kutta_dopri5< state_type > dopri5_type; |
| typedef controlled_runge_kutta< dopri5_type > controlled_dopri5_type; |
| typedef dense_output_runge_kutta< controlled_dopri5_type > dense_output_dopri5_type; |
| |
| dense_output_dopri5_type dopri5 = make_dense_output( 1E-9 , 1E-9 , dopri5_type() ); |
| |
| x[0] = 2.0 / sqrt(3.0); |
| |
| out.open( "bs4.dat" ); |
| out.precision(16); |
| integrate_adaptive( dopri5 , rhs , x , t , t_end , dt , write_out ); |
| out.close(); |
| |
| } |