123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147 |
- /*
- 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;
- }
|