1 /*
2 libs/numeric/odeint/examples/stochastic_euler.hpp
3
4 Copyright 2012 Karsten Ahnert
5 Copyright 2012 Mario Mulansky
6
7 Stochastic euler stepper example and Ornstein-Uhlenbeck process
8
9 Distributed under the Boost Software License, Version 1.0.
10 (See accompanying file LICENSE_1_0.txt or
11 copy at http://www.boost.org/LICENSE_1_0.txt)
12 */
13
14
15 #include <vector>
16 #include <iostream>
17 #include <boost/random.hpp>
18 #include <boost/array.hpp>
19
20 #include <boost/numeric/odeint.hpp>
21
22
23 /*
24 //[ stochastic_euler_class_definition
25 template< size_t N > class stochastic_euler
26 {
27 public:
28
29 typedef boost::array< double , N > state_type;
30 typedef boost::array< double , N > deriv_type;
31 typedef double value_type;
32 typedef double time_type;
33 typedef unsigned short order_type;
34 typedef boost::numeric::odeint::stepper_tag stepper_category;
35
36 static order_type order( void ) { return 1; }
37
38 // ...
39 };
40 //]
41 */
42
43
44 /*
45 //[ stochastic_euler_do_step
46 template< size_t N > class stochastic_euler
47 {
48 public:
49
50 // ...
51
52 template< class System >
53 void do_step( System system , state_type &x , time_type t , time_type dt ) const
54 {
55 deriv_type det , stoch ;
56 system.first( x , det );
57 system.second( x , stoch );
58 for( size_t i=0 ; i<x.size() ; ++i )
59 x[i] += dt * det[i] + sqrt( dt ) * stoch[i];
60 }
61 };
62 //]
63 */
64
65
66
67
68 //[ stochastic_euler_class
69 template< size_t N >
70 class stochastic_euler
71 {
72 public:
73
74 typedef boost::array< double , N > state_type;
75 typedef boost::array< double , N > deriv_type;
76 typedef double value_type;
77 typedef double time_type;
78 typedef unsigned short order_type;
79
80 typedef boost::numeric::odeint::stepper_tag stepper_category;
81
order(void)82 static order_type order( void ) { return 1; }
83
84 template< class System >
do_step(System system,state_type & x,time_type t,time_type dt) const85 void do_step( System system , state_type &x , time_type t , time_type dt ) const
86 {
87 deriv_type det , stoch ;
88 system.first( x , det );
89 system.second( x , stoch );
90 for( size_t i=0 ; i<x.size() ; ++i )
91 x[i] += dt * det[i] + sqrt( dt ) * stoch[i];
92 }
93 };
94 //]
95
96
97
98 //[ stochastic_euler_ornstein_uhlenbeck_def
99 const static size_t N = 1;
100 typedef boost::array< double , N > state_type;
101
102 struct ornstein_det
103 {
operator ()ornstein_det104 void operator()( const state_type &x , state_type &dxdt ) const
105 {
106 dxdt[0] = -x[0];
107 }
108 };
109
110 struct ornstein_stoch
111 {
112 boost::mt19937 &m_rng;
113 boost::normal_distribution<> m_dist;
114
ornstein_stochornstein_stoch115 ornstein_stoch( boost::mt19937 &rng , double sigma ) : m_rng( rng ) , m_dist( 0.0 , sigma ) { }
116
operator ()ornstein_stoch117 void operator()( const state_type &x , state_type &dxdt )
118 {
119 dxdt[0] = m_dist( m_rng );
120 }
121 };
122 //]
123
124 struct streaming_observer
125 {
126 template< class State >
operator ()streaming_observer127 void operator()( const State &x , double t ) const
128 {
129 std::cout << t << "\t" << x[0] << "\n";
130 }
131 };
132
133
main(int argc,char ** argv)134 int main( int argc , char **argv )
135 {
136 using namespace std;
137 using namespace boost::numeric::odeint;
138
139 //[ ornstein_uhlenbeck_main
140 boost::mt19937 rng;
141 double dt = 0.1;
142 state_type x = {{ 1.0 }};
143 integrate_const( stochastic_euler< N >() , make_pair( ornstein_det() , ornstein_stoch( rng , 1.0 ) ),
144 x , 0.0 , 10.0 , dt , streaming_observer() );
145 //]
146 return 0;
147 }
148