| /* |
| * stepper_details.cpp |
| * |
| * This example demonstrates some details about the steppers in odeint. |
| * |
| * |
| * Copyright 2011-2012 Karsten Ahnert |
| * Copyright 2012 Mario Mulansky |
| * Copyright 2013 Pascal Germroth |
| * |
| * 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 <boost/array.hpp> |
| #include <boost/bind.hpp> |
| #include <boost/numeric/odeint.hpp> |
| |
| using namespace std; |
| using namespace boost::numeric::odeint; |
| |
| const size_t N = 3; |
| |
| typedef boost::array< double , N > state_type; |
| |
| //[ system_function_structure |
| void sys( const state_type & /*x*/ , state_type & /*dxdt*/ , const double /*t*/ ) |
| { |
| // ... |
| } |
| //] |
| |
| void sys1( const state_type &/*x*/ , state_type &/*dxdt*/ , const double /*t*/ ) |
| { |
| } |
| |
| void sys2( const state_type &/*x*/ , state_type &/*dxdt*/ , const double /*t*/ ) |
| { |
| } |
| |
| |
| //[ symplectic_stepper_detail_system_function |
| typedef boost::array< double , 1 > vector_type; |
| |
| |
| struct harm_osc_f1 |
| { |
| void operator()( const vector_type &p , vector_type &dqdt ) |
| { |
| dqdt[0] = p[0]; |
| } |
| }; |
| |
| struct harm_osc_f2 |
| { |
| void operator()( const vector_type &q , vector_type &dpdt ) |
| { |
| dpdt[0] = -q[0]; |
| } |
| }; |
| //] |
| |
| //[ symplectic_stepper_detail_system_class |
| struct harm_osc |
| { |
| void f1( const vector_type &p , vector_type &dqdt ) const |
| { |
| dqdt[0] = p[0]; |
| } |
| |
| void f2( const vector_type &q , vector_type &dpdt ) const |
| { |
| dpdt[0] = -q[0]; |
| } |
| }; |
| //] |
| |
| int main( int argc , char **argv ) |
| { |
| using namespace std; |
| |
| // Explicit stepper example |
| { |
| double t( 0.0 ) , dt( 0.1 ); |
| state_type in , out , dxdtin , inout; |
| //[ explicit_stepper_detail_example |
| runge_kutta4< state_type > rk; |
| rk.do_step( sys1 , inout , t , dt ); // In-place transformation of inout |
| rk.do_step( sys2 , inout , t , dt ); // call with different system: Ok |
| rk.do_step( sys1 , in , t , out , dt ); // Out-of-place transformation |
| rk.do_step( sys1 , inout , dxdtin , t , dt ); // In-place tranformation of inout |
| rk.do_step( sys1 , in , dxdtin , t , out , dt ); // Out-of-place transformation |
| //] |
| } |
| |
| |
| |
| // FSAL stepper example |
| { |
| double t( 0.0 ) , dt( 0.1 ); |
| state_type in , in2 , in3 , out , dxdtin , dxdtout , inout , dxdtinout; |
| //[ fsal_stepper_detail_example |
| runge_kutta_dopri5< state_type > rk; |
| rk.do_step( sys1 , in , t , out , dt ); |
| rk.do_step( sys2 , in , t , out , dt ); // DONT do this, sys1 is assumed |
| |
| rk.do_step( sys2 , in2 , t , out , dt ); |
| rk.do_step( sys2 , in3 , t , out , dt ); // DONT do this, in2 is assumed |
| |
| rk.do_step( sys1 , inout , dxdtinout , t , dt ); |
| rk.do_step( sys2 , inout , dxdtinout , t , dt ); // Ok, internal derivative is not used, dxdtinout is updated |
| |
| rk.do_step( sys1 , in , dxdtin , t , out , dxdtout , dt ); |
| rk.do_step( sys2 , in , dxdtin , t , out , dxdtout , dt ); // Ok, internal derivative is not used |
| //] |
| } |
| |
| |
| // Symplectic harmonic oscillator example |
| { |
| double t( 0.0 ) , dt( 0.1 ); |
| //[ symplectic_stepper_detail_example |
| pair< vector_type , vector_type > x; |
| x.first[0] = 1.0; x.second[0] = 0.0; |
| symplectic_rkn_sb3a_mclachlan< vector_type > rkn; |
| rkn.do_step( make_pair( harm_osc_f1() , harm_osc_f2() ) , x , t , dt ); |
| //] |
| |
| //[ symplectic_stepper_detail_system_class_example |
| harm_osc h; |
| rkn.do_step( make_pair( boost::bind( &harm_osc::f1 , h , _1 , _2 ) , boost::bind( &harm_osc::f2 , h , _1 , _2 ) ) , |
| x , t , dt ); |
| //] |
| } |
| |
| // Simplified harmonic oscillator example |
| { |
| double t = 0.0, dt = 0.1; |
| //[ simplified_symplectic_stepper_example |
| pair< vector_type , vector_type > x; |
| x.first[0] = 1.0; x.second[0] = 0.0; |
| symplectic_rkn_sb3a_mclachlan< vector_type > rkn; |
| rkn.do_step( harm_osc_f1() , x , t , dt ); |
| //] |
| |
| vector_type q = {{ 1.0 }} , p = {{ 0.0 }}; |
| //[ symplectic_stepper_detail_ref_usage |
| rkn.do_step( harm_osc_f1() , make_pair( boost::ref( q ) , boost::ref( p ) ) , t , dt ); |
| rkn.do_step( harm_osc_f1() , q , p , t , dt ); |
| rkn.do_step( make_pair( harm_osc_f1() , harm_osc_f2() ) , q , p , t , dt ); |
| //] |
| } |
| |
| // adams_bashforth_moulton stepper example |
| { |
| double t = 0.0 , dt = 0.1; |
| state_type inout; |
| //[ multistep_detail_example |
| adams_bashforth_moulton< 5 , state_type > abm; |
| abm.initialize( sys , inout , t , dt ); |
| abm.do_step( sys , inout , t , dt ); |
| //] |
| |
| //[ multistep_detail_own_stepper_initialization |
| abm.initialize( runge_kutta_fehlberg78< state_type >() , sys , inout , t , dt ); |
| //] |
| } |
| |
| |
| |
| // dense output stepper examples |
| { |
| double t = 0.0 , dt = 0.1; |
| state_type in; |
| //[ dense_output_detail_example |
| dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > > dense; |
| dense.initialize( in , t , dt ); |
| pair< double , double > times = dense.do_step( sys ); |
| (void)times; |
| //] |
| |
| state_type inout; |
| double t_start = 0.0 , t_end = 1.0; |
| //[ dense_output_detail_generation1 |
| typedef boost::numeric::odeint::result_of::make_dense_output< |
| runge_kutta_dopri5< state_type > >::type dense_stepper_type; |
| dense_stepper_type dense2 = make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() ); |
| (void)dense2; |
| //] |
| |
| //[ dense_output_detail_generation2 |
| integrate_const( make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() ) , sys , inout , t_start , t_end , dt ); |
| //] |
| } |
| |
| |
| return 0; |
| } |