From 83c55a46326ab268cca3de8a3bf1c482de8ffd78 Mon Sep 17 00:00:00 2001
From: Karsten Ahnert
where ξ is Gaussian white noise with zero mean and
- a standard deviation σ. f(x) is
- said to be the deterministic part while g(x) ξ is the
- noisy part. In case g(x) is independent of x
+ a standard deviation σ(t). f(x)
+ is said to be the deterministic part while g(x) ξ is
+ the noisy part. In case g(x) is independent of x
the SDE is said to have additive noise. It is not possible to solve SDE
with the classical solvers for ODEs since the noisy part of the SDE has
to be scaled differently then the deterministic part with respect to the
@@ -1964,8 +1964,7 @@
deterministic and one for the stochastic part. The first element of the
pair simply computes the deterministic part while the second the stochastic
one. Then, the second part also needs to calculate the random numbers in
- order to simulate the stochastic process. We can now implement the do_step
- method
+ order to simulate the stochastic process. We can now implement the do_step
method
@@ -2009,7 +2008,65 @@
- examples, with integrate function + Now, lets look how we use the new stepper. A nice example is the Ornstein-Uhlenbeck + process. It consists of a simple Brownian motion overlapped with an relaxation + process. Its SDE reads +
++ dx/dt = - x + ξ +
++ where ξ is Gaussian white noise with standard deviation σ. + Implementing the Ornstein-Uhlenbeck process is quite simple. We need two + functions or funtors - one for the deterministic and one for the stochastic + part: +
++
+const static size_t N = 1; +typedef std::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 +{ + std::mt19937 m_rng; + std::normal_distribution<> m_dist; + + ornstein_stoch( double sigma ) : m_rng() , m_dist( 0.0 , sigma ) { } + + void operator()( const state_type &x , state_type &dxdt ) + { + dxdt[0] = m_dist( m_rng ); + } +}; ++
+
+
+ In the stochastic part we have used the Mersenne twister for the random
+ number generation and a Gaussian white noise generator normal_distribution
+ with standard deviation σ. Now, we can use the stochastic
+ Euler stepper with the integrate functions:
+
+
+double dt = 0.1; +state_type x = {{ 1.0 }}; +integrate_const( stochastic_euler< N >() , make_pair( ornstein_det() , ornstein_stoch( 1.0 ) ) , + x , 0.0 , 10.0 , dt , streaming_observer() ); ++
+
+
+ Note, how we have used the make_pair
+ function for the generation of the system function.
+ Most steppers in odeint also accept the state give as a range. A range is + sequence of values modelled by a range concept, see Boost.Range + for an overview over existing concepts and examples of ranges. +
++ state type mus not necessary be used in the steppers +
++ example +
++ table which steppers support boost::range and with which algebra +
+