| /* |
| libs/numeric/odeint/examples/stochastic_euler.hpp |
| |
| Copyright 2012 Karsten Ahnert |
| Copyright 2012 Mario Mulansky |
| |
| Stochastic euler stepper example and Ornstein-Uhlenbeck process |
| |
| 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 <vector> |
| #include <iostream> |
| #include <boost/random.hpp> |
| #include <boost/array.hpp> |
| |
| #include <boost/numeric/odeint.hpp> |
| |
| |
| /* |
| //[ stochastic_euler_class_definition |
| template< size_t N > class stochastic_euler |
| { |
| public: |
| |
| typedef boost::array< double , N > state_type; |
| typedef boost::array< double , N > deriv_type; |
| typedef double value_type; |
| typedef double time_type; |
| typedef unsigned short order_type; |
| typedef boost::numeric::odeint::stepper_tag stepper_category; |
| |
| static order_type order( void ) { return 1; } |
| |
| // ... |
| }; |
| //] |
| */ |
| |
| |
| /* |
| //[ stochastic_euler_do_step |
| template< size_t N > class stochastic_euler |
| { |
| public: |
| |
| // ... |
| |
| template< class System > |
| void do_step( System system , state_type &x , time_type t , time_type dt ) const |
| { |
| deriv_type det , stoch ; |
| system.first( x , det ); |
| system.second( x , stoch ); |
| for( size_t i=0 ; i<x.size() ; ++i ) |
| x[i] += dt * det[i] + sqrt( dt ) * stoch[i]; |
| } |
| }; |
| //] |
| */ |
| |
| |
| |
| |
| //[ stochastic_euler_class |
| template< size_t N > |
| class stochastic_euler |
| { |
| public: |
| |
| typedef boost::array< double , N > state_type; |
| typedef boost::array< double , N > deriv_type; |
| typedef double value_type; |
| typedef double time_type; |
| typedef unsigned short order_type; |
| |
| typedef boost::numeric::odeint::stepper_tag stepper_category; |
| |
| static order_type order( void ) { return 1; } |
| |
| template< class System > |
| void do_step( System system , state_type &x , time_type t , time_type dt ) const |
| { |
| deriv_type det , stoch ; |
| system.first( x , det ); |
| system.second( x , stoch ); |
| for( size_t i=0 ; i<x.size() ; ++i ) |
| x[i] += dt * det[i] + sqrt( dt ) * stoch[i]; |
| } |
| }; |
| //] |
| |
| |
| |
| //[ stochastic_euler_ornstein_uhlenbeck_def |
| const static size_t N = 1; |
| typedef boost::array< double , N > state_type; |
| |
| struct ornstein_det |
| { |
| void operator()( const state_type &x , state_type &dxdt ) const |
| { |
| dxdt[0] = -x[0]; |
| } |
| }; |
| |
| struct ornstein_stoch |
| { |
| boost::mt19937 &m_rng; |
| boost::normal_distribution<> m_dist; |
| |
| ornstein_stoch( boost::mt19937 &rng , double sigma ) : m_rng( rng ) , m_dist( 0.0 , sigma ) { } |
| |
| void operator()( const state_type &x , state_type &dxdt ) |
| { |
| dxdt[0] = m_dist( m_rng ); |
| } |
| }; |
| //] |
| |
| struct streaming_observer |
| { |
| template< class State > |
| void operator()( const State &x , double t ) const |
| { |
| std::cout << t << "\t" << x[0] << "\n"; |
| } |
| }; |
| |
| |
| int main( int argc , char **argv ) |
| { |
| using namespace std; |
| using namespace boost::numeric::odeint; |
| |
| //[ ornstein_uhlenbeck_main |
| boost::mt19937 rng; |
| double dt = 0.1; |
| state_type x = {{ 1.0 }}; |
| integrate_const( stochastic_euler< N >() , make_pair( ornstein_det() , ornstein_stoch( rng , 1.0 ) ), |
| x , 0.0 , 10.0 , dt , streaming_observer() ); |
| //] |
| return 0; |
| } |