diff --git a/doc/boost_sandbox_numeric_odeint/concepts/controlled_stepper.html b/doc/boost_sandbox_numeric_odeint/concepts/controlled_stepper.html index 37af938a..f5e3cc46 100644 --- a/doc/boost_sandbox_numeric_odeint/concepts/controlled_stepper.html +++ b/doc/boost_sandbox_numeric_odeint/concepts/controlled_stepper.html @@ -25,9 +25,9 @@ to be used within integrate functions.

-
- - Description +
+ + Description

A controlled stepper following this Controlled Stepper concept provides the @@ -36,9 +36,9 @@ with a given step-size dt. Depending on an error estimate of the solution the step might be rejected and a smaller step-size is suggested.

-
- - Notation +
+ + Notation

@@ -77,9 +77,9 @@

-
- - Valid +
+ + Valid Expressions
@@ -178,9 +178,9 @@
-
- - Models +
+ + Models
-
- - Description +
+ + Description

This concept describes how to define a ODE that can be solved by an implicit @@ -33,9 +33,9 @@ therefore, the state_type implicit routines is ublas::vector and the matrix is defined as ublas::matrix.

-
- - Notation +
+ + Notation

@@ -70,9 +70,9 @@

-
- - Valid +
+ + Valid Expressions
diff --git a/doc/boost_sandbox_numeric_odeint/concepts/observer.html b/doc/boost_sandbox_numeric_odeint/concepts/observer.html index 1a8ba9f4..c58f79c1 100644 --- a/doc/boost_sandbox_numeric_odeint/concepts/observer.html +++ b/doc/boost_sandbox_numeric_odeint/concepts/observer.html @@ -19,9 +19,9 @@ -
- - Description +
+ + Description

The Observer concept is used in integrate @@ -32,9 +32,9 @@ or at user defined time points depending on which integrate function is used.

-
- - Notation +
+ + Notation

@@ -65,9 +65,9 @@

-
- - Valid +
+ + Valid Expressions
diff --git a/doc/boost_sandbox_numeric_odeint/concepts/simple_symplectic_system.html b/doc/boost_sandbox_numeric_odeint/concepts/simple_symplectic_system.html index 1e830be6..ca0ef162 100644 --- a/doc/boost_sandbox_numeric_odeint/concepts/simple_symplectic_system.html +++ b/doc/boost_sandbox_numeric_odeint/concepts/simple_symplectic_system.html @@ -20,9 +20,9 @@ Simple Symplectic System -
- - Description +
+ + Description

In most Hamiltonian systems the kinetic term is a quadratic term in the momentum @@ -44,9 +44,9 @@ equation for p' has to be provided to the symplectic stepper. We call this concept Simple_Symplectic_System

-
- - Notation +
+ + Notation

@@ -77,9 +77,9 @@

-
- - Valid +
+ + Valid Expressions
diff --git a/doc/boost_sandbox_numeric_odeint/concepts/state_algebra_operations.html b/doc/boost_sandbox_numeric_odeint/concepts/state_algebra_operations.html index 52379952..3e61e5c6 100644 --- a/doc/boost_sandbox_numeric_odeint/concepts/state_algebra_operations.html +++ b/doc/boost_sandbox_numeric_odeint/concepts/state_algebra_operations.html @@ -40,9 +40,9 @@ and no algebra/operations are used.

-
- - Description +
+ + Description

The State, Algebra and Operations @@ -82,9 +82,9 @@

Operations

-
- - Notation +
+ + Notation

@@ -142,9 +142,9 @@

-
- - Valid +
+ + Valid Expressions
@@ -414,9 +414,9 @@ -
- - Notation +
+ + Notation

@@ -452,9 +452,9 @@

-
- - Valid +
+ + Valid Expressions
diff --git a/doc/boost_sandbox_numeric_odeint/concepts/state_wrapper.html b/doc/boost_sandbox_numeric_odeint/concepts/state_wrapper.html index 9733f965..587376e0 100644 --- a/doc/boost_sandbox_numeric_odeint/concepts/state_wrapper.html +++ b/doc/boost_sandbox_numeric_odeint/concepts/state_wrapper.html @@ -20,9 +20,9 @@ State Wrapper -
- - Description +
+ + Description

The State Wrapper @@ -30,9 +30,9 @@ intermediate results within the stepper's do_step methods.

-
- - Notation +
+ + Notation

@@ -56,9 +56,9 @@

-
- - Valid +
+ + Valid Expressions
diff --git a/doc/boost_sandbox_numeric_odeint/concepts/stepper.html b/doc/boost_sandbox_numeric_odeint/concepts/stepper.html index 3279fc2b..c4633fb5 100644 --- a/doc/boost_sandbox_numeric_odeint/concepts/stepper.html +++ b/doc/boost_sandbox_numeric_odeint/concepts/stepper.html @@ -24,9 +24,9 @@ be used within the integrate functions.

-
- - Description +
+ + Description

The basic stepper concept. A basic stepper following this Stepper concept @@ -40,9 +40,9 @@ Symplectic System or Implicit System. Note that all error steppers are also basic steppers.

-
- - Refinement +
+ + Refinement of
    @@ -53,9 +53,9 @@ CopyConstructable
-
- - Notation +
+ + Notation

@@ -95,9 +95,9 @@

-
- - Valid +
+ + Valid Expressions
@@ -187,9 +187,9 @@
-
- - Models +
+ + Models
-
- - Description +
+ + Description

This concept describes how to define a symplectic system written with generalized @@ -58,9 +58,9 @@ std::pair of functions or a functors. Symplectic systems are used in symplectic steppers like symplectic_rkn_sb3a_mclachlan.

-
- - Notation +
+ + Notation

@@ -107,9 +107,9 @@

-
- - Valid +
+ + Valid expressions
diff --git a/doc/boost_sandbox_numeric_odeint/concepts/system.html b/doc/boost_sandbox_numeric_odeint/concepts/system.html index 8c66523e..d3e50519 100644 --- a/doc/boost_sandbox_numeric_odeint/concepts/system.html +++ b/doc/boost_sandbox_numeric_odeint/concepts/system.html @@ -19,9 +19,9 @@ -
- - Description +
+ + Description

The System concept models the algorithmic implementation of the rhs. of the @@ -34,9 +34,9 @@ System and Implicit System.

-
- - Notation +
+ + Notation

@@ -75,9 +75,9 @@

-
- - Valid +
+ + Valid expressions
diff --git a/doc/boost_sandbox_numeric_odeint/getting_started/short_example.html b/doc/boost_sandbox_numeric_odeint/getting_started/short_example.html index a939f168..e92f2639 100644 --- a/doc/boost_sandbox_numeric_odeint/getting_started/short_example.html +++ b/doc/boost_sandbox_numeric_odeint/getting_started/short_example.html @@ -30,7 +30,6 @@ = (p,-q+γ p):

-

/* The type of container used to hold the state vector */
 typedef std::vector< double > state_type;
@@ -60,7 +59,6 @@
         start:
       

-

state_type x(2);
 x[0] = 1.0; // start at x=1.0, p=0.0
@@ -75,7 +73,6 @@
         stepper (5th order) and uses adaptive stepsize.
       

-

size_t steps = integrate( harmonic_oscillator ,
         x , 0.0 , 10.0 , 0.1 );
@@ -96,7 +93,6 @@
         rhs must then be implemented as a functor having defined the ()-operator:
       

-

/* The rhs of x' = f(x) defined as a class */
 class harm_osc {
@@ -119,7 +115,6 @@
         which can be used via
       

-

harm_osc ho(0.15);
 steps = integrate( ho ,
@@ -134,7 +129,6 @@
         to provide a reasonable observer. An example is
       

-

struct push_back_state_and_time
 {
@@ -158,7 +152,6 @@
         pass this container to the integration function:
       

-

vector<state_type> x_vec;
 vector<double> times;
diff --git a/doc/boost_sandbox_numeric_odeint/getting_started/usage__compilation__headers.html b/doc/boost_sandbox_numeric_odeint/getting_started/usage__compilation__headers.html
index 04b5f298..81c1cb6e 100644
--- a/doc/boost_sandbox_numeric_odeint/getting_started/usage__compilation__headers.html
+++ b/doc/boost_sandbox_numeric_odeint/getting_started/usage__compilation__headers.html
@@ -25,13 +25,12 @@
         library. It can be include by
       

-

#include <boost/numeric/odeint.hpp>
 

which includes all headers in the library. All functions and classes from - odeint live in the namespace + odeint live in the namespace

using namespace boost::numeric::odeint;
 
diff --git a/doc/boost_sandbox_numeric_odeint/odeint_in_detail/integrate_functions.html b/doc/boost_sandbox_numeric_odeint/odeint_in_detail/integrate_functions.html index c8c2e888..a70d7a85 100644 --- a/doc/boost_sandbox_numeric_odeint/odeint_in_detail/integrate_functions.html +++ b/doc/boost_sandbox_numeric_odeint/odeint_in_detail/integrate_functions.html @@ -35,9 +35,9 @@ Output Stepper. Depending on the abilities of the stepper, the integrate functions make use of step-size control or dense output.

-
- - Equidistant +
+ + Equidistant observer calls

@@ -109,9 +109,9 @@ dt. -

- - Observer +
+ + Observer calls at each step

@@ -179,9 +179,9 @@ Stepper above. No dense output is used. -

- - Observer +
+ + Observer calls at given time points

@@ -254,9 +254,9 @@ x(t) at the time points from the sequence. -

- - Integrate +
+ + Integrate a given number of steps

@@ -311,9 +311,9 @@ This function returns the end time t0 + n*dt.

-
- - Convenience +
+ + Convenience integrate function

diff --git a/doc/boost_sandbox_numeric_odeint/odeint_in_detail/state_types__algebras_and_operations.html b/doc/boost_sandbox_numeric_odeint/odeint_in_detail/state_types__algebras_and_operations.html index de6957c2..c3871cc3 100644 --- a/doc/boost_sandbox_numeric_odeint/odeint_in_detail/state_types__algebras_and_operations.html +++ b/doc/boost_sandbox_numeric_odeint/odeint_in_detail/state_types__algebras_and_operations.html @@ -399,7 +399,6 @@ as resizeable as it is a dynamically sized array:

-

class my_vector : public std::vector< double >
 {
@@ -534,7 +533,6 @@
               specializations:
             

-

typedef std::list< double > state_type;
 
@@ -547,7 +545,7 @@
     const static bool value = type::value;
 };
 
-template< > 
+template< >
 struct same_size_impl< state_type , state_type >
 { // define how to check size
     static bool same_size( const state_type &v1 ,
@@ -567,7 +565,7 @@
     }
 };
 
-} } } 
+} } }
 

@@ -627,7 +625,6 @@ instances of state types:

-

template<>
 struct state_wrapper< gsl_vector* >
@@ -664,7 +661,7 @@
               This state_wrapper
               specialization tells odeint how gsl_vectors are created, copied and
               destroyed. Next we need resizing, this is required because gsl_vectors
-              are dynamically sized objects: 
+              are dynamically sized objects:
 

template<>
 struct is_resizeable< gsl_vector* >
@@ -699,13 +696,12 @@
               support iterators, so we first implement a gsl iterator:
             

-

/*
  * defines an iterator for gsl_vector
  */
-class gsl_vector_iterator 
-      : public boost::iterator_facade< gsl_vector_iterator , double , 
+class gsl_vector_iterator
+      : public boost::iterator_facade< gsl_vector_iterator , double ,
                                        boost::random_access_traversal_tag >
 {
 public :
@@ -734,7 +730,7 @@
               A similar class exists for the const
               version of the iterator. Then we have a function returning the end
               iterator (similarily for const
-              again): 
+              again):
 

gsl_vector_iterator end_iterator( gsl_vector *x )
 {
@@ -747,7 +743,7 @@
             

Finally, the bindings for Boost.Range - are added: + are added:

// template<>
 inline gsl_vector_iterator range_begin( gsl_vector *x )
@@ -1061,7 +1057,6 @@
               function so all that has to be done in this case is to declare resizability:
             

-

typedef boost::numeric::ublas::vector< double > state_type;
 
@@ -1085,7 +1080,6 @@
               code shows the corresponding definitions:
             

-

int main()
 {
@@ -1118,9 +1112,8 @@
               type looks as follows:
             

-

-
class point3D :     
+
class point3D :
     boost::additive1< point3D ,
     boost::additive2< point3D , double ,
     boost::multiplicative2< point3D , double > > >
@@ -1132,7 +1125,7 @@
     point3D()
         : x( 0.0 ) , y( 0.0 ) , z( 0.0 )
     { }
-    
+
     point3D( const double val )
         : x( val ) , y( val ) , z( val )
     { }
@@ -1173,7 +1166,6 @@
               type:
             

-

// only required for steppers with error control
 point3D operator/( const point3D &p1 , const point3D &p2 )
@@ -1193,7 +1185,6 @@
               implementing a reduction over the state type:
             

-

namespace boost { namespace numeric { namespace odeint {
 // specialization of vector_space_reduce, only required for steppers with error control
@@ -1209,7 +1200,7 @@
         return init;
     }
 };
-} } } 
+} } }
 

@@ -1221,7 +1212,6 @@ easily perform the integration on a lorenz system by using the vector_space_algebra again:

-

const double sigma = 10.0;
 const double R = 28.0;
@@ -1242,9 +1232,9 @@
     point3D x( 10.0 , 5.0 , 5.0 );
     //typedef runge_kutta4< point3D , double , point3D , 
     //                      double , vector_space_algebra > stepper;
-    typedef runge_kutta_dopri5< point3D , double , point3D , 
+    typedef runge_kutta_dopri5< point3D , double , point3D ,
                                 double , vector_space_algebra > stepper;
-    integrate_adaptive( make_controlled<stepper>( 1E-8 , 1E-8 ) , lorenz , x , 
+    integrate_adaptive( make_controlled<stepper>( 1E-8 , 1E-8 ) , lorenz , x ,
                         0.0 , 10.0 , 0.1 );
     std::cout << x << std::endl;
 }
diff --git a/doc/boost_sandbox_numeric_odeint/odeint_in_detail/steppers.html b/doc/boost_sandbox_numeric_odeint/odeint_in_detail/steppers.html
index 9f92349a..951f5431 100644
--- a/doc/boost_sandbox_numeric_odeint/odeint_in_detail/steppers.html
+++ b/doc/boost_sandbox_numeric_odeint/odeint_in_detail/steppers.html
@@ -132,7 +132,6 @@
           is the Runge Kutta stepper of fourth order:
         

-

runge_kutta4< state_type > rk;
 rk.do_step( sys1 , inout , t , dt );               // In-place transformation of inout
@@ -177,7 +176,6 @@
           how the FSAL steppers can be used is
         

-

runge_kutta_dopri5< state_type > rk;
 rk.do_step( sys1 , in , t , out , dt );
@@ -243,7 +241,6 @@
           again simple C-functions of functors. An example is the harmonic oscillator:
         

-

typedef boost::array< double , 1 > vector_type;
 
@@ -272,7 +269,6 @@
           example for the harmonic oscillator is now:
         

-

pair< vector_type , vector_type > x;
 x.first[0] = 1.0; x.second[0] = 0.0;
@@ -286,7 +282,6 @@
           two public method:
         

-

struct harm_osc
 {
@@ -304,7 +299,6 @@
 

-

harm_osc h;
 rkn.do_step( make_pair( boost::bind( &harm_osc::f1 , h , _1 , _2 ) , boost::bind( &harm_osc::f2 , h , _1 , _2 ) ) ,
@@ -324,7 +318,6 @@
           can be written as
         

-

pair< vector_type , vector_type > x;
 x.first[0] = 1.0; x.second[0] = 0.0;
@@ -348,7 +341,6 @@
           takes q and p without combining them into a pair:
         

-

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 );
@@ -414,7 +406,6 @@
           steps is instantiated and initialized;
         

-

adams_bashforth_moulton< 5 , state_type > abm;
 abm.initialize( sys , inout , t , dt );
@@ -432,7 +423,6 @@
           state of the Adams-Bashforth-Stepper:
         

-

abm.initialize( runge_kutta_fehlberg78< state_type >() , sys , inout , t , dt );
 
@@ -695,7 +685,6 @@ An example is

-

dense_output_runge_kutta< controlled_runge_kutta< runge_kutta_dopri5< state_type > > > dense;
 dense.initialize( in , t , dt );
@@ -721,7 +710,6 @@
           exist:
         

-

boost::numeric::odeint::result_of::make_dense_output< runge_kutta_dopri5< state_type > >::type dense2 = make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() );
 
@@ -739,7 +727,6 @@ for easy use with the integrate functions:

-

integrate_const( make_dense_output( 1.0e-6 , 1.0e-6 , runge_kutta_dopri5< state_type >() ) , sys , inout , t_start , t_end , dt );
 
@@ -808,7 +795,6 @@ 4 internal derivatives of the solution at times (t-dt,t-2dt,t-3dt,t-4dt).

-

adams_bashforth_moulton< 4 , state_type > stepper;
 stepper.do_step( sys , x , t , dt );   // make one step with the classical Runge-Kutta stepper and initialize the first internal state at time t-4dt
@@ -1939,39 +1925,89 @@
           x(t+Δ t) = x(t) + Δ t f(x(t)) + Δ t1/2 g(x) ξ(t)
         

- where ξ(t) is a independent and identically distributed normal random variable. + where ξ(t) is an independent normal distributed random variable.

- Now we will implement this method. We will call the stepper stochastic_euler. Its class definition - looks like + Now we will implement this method. We will call the stepper stochastic_euler. It models the __stepper_concept. + The class definition looks like

-

-
class stochastic_euler
+
template< size_t N > class stochastic_euler
 {
 public:
 
-    typedef std::vector< double > state_type;
-    typedef std::vector< double > deriv_type;
+    typedef std::array< double , N > state_type;
+    typedef std::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; }
+
+    // ...
 };
 

- System function + The types are needed in order to fulfill the stepper concept. As internal + state and deriv type we use simple arrays in the stochastic euler, they + are needed for the temporaries. The stepper has the order one which is + returned from the order() function.

- steppers body + The system functions needs to calculate the deterministic and the stochastic + part of our stochastic differential equation. So it might be suitable that + the system function itself is a pair of functions, one for computing the + 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

- do_step method +

+
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];
+    }
+};
+
+

+

+ This is all. It is quite simple and the stochastic euler stepper implement + here is quite general. Of course it can be enhanced, for example +

+
    +
  • + use of operations and algebras as well as the resizing mechanism for + maximal flexibility and portability +
  • +
  • + use of boost::ref for the system functions +
  • +
  • + use of boost::range for the state type in the + do_step method +
  • +
  • + ... +
  • +

examples, with integrate function

diff --git a/doc/boost_sandbox_numeric_odeint/tutorial/chaotic_systems_and_lyapunov_exponents.html b/doc/boost_sandbox_numeric_odeint/tutorial/chaotic_systems_and_lyapunov_exponents.html index f576b099..566791b9 100644 --- a/doc/boost_sandbox_numeric_odeint/tutorial/chaotic_systems_and_lyapunov_exponents.html +++ b/doc/boost_sandbox_numeric_odeint/tutorial/chaotic_systems_and_lyapunov_exponents.html @@ -61,7 +61,6 @@ The implementation of the Lorenz system is

-

const double sigma = 10.0;
 const double R = 28.0;
@@ -84,7 +83,6 @@
         the perturbation could look like:
       

-

const size_t n = 3;
 const size_t num_of_lyap = 3;
@@ -116,7 +114,6 @@
         inside a functor with templatized arguments:
       

-

struct lorenz
 {
@@ -138,7 +135,7 @@
 

This works fine and lorenz_with_lyap - can be used for example via + can be used for example via

state_type x;
 // initialize x
@@ -154,7 +151,6 @@
         current state lies on the attractor, hence it would look like
       

-

state_type x;
 // initialize x
@@ -170,7 +166,6 @@
         by pointer to the state and an appropriate length, something similar to
       

-

void lorenz( double* x , double *dxdt , double t, void* params )
 {
@@ -188,7 +183,6 @@
         the system has to by changed:
       

-

struct lorenz
 {
@@ -211,7 +205,6 @@
         with an suited range:
       

-

// perform 10000 transient steps
 integrate_n_steps( rk4 , lorenz() , std::make_pair( x.begin() , x.begin() + n ) , 0.0 , dt , 10000 );
@@ -242,7 +235,6 @@
           
 
 

-

fill( x.begin()+n , x.end() , 0.0 );
 for( size_t i=0 ; i<num_of_lyap ; ++i ) x[n+n*i+i] = 1.0;
diff --git a/doc/boost_sandbox_numeric_odeint/tutorial/harmonic_oscillator.html b/doc/boost_sandbox_numeric_odeint/tutorial/harmonic_oscillator.html
index 0836da63..5332e1bd 100644
--- a/doc/boost_sandbox_numeric_odeint/tutorial/harmonic_oscillator.html
+++ b/doc/boost_sandbox_numeric_odeint/tutorial/harmonic_oscillator.html
@@ -52,7 +52,6 @@
           to just define a function, e.g:
         

-

/* The type of container used to hold the state vector */
 typedef std::vector< double > state_type;
@@ -79,7 +78,6 @@
           parameter structure as above:
         

-

/* The rhs of x' = f(x) defined as a class */
 class harm_osc {
@@ -1138,7 +1136,6 @@
           ) function from odeint:
         

-

runge_kutta4< state_type > stepper;
 integrate_const( stepper , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
@@ -1155,7 +1152,6 @@
           method which can used directly. So, you write down the above example as
         

-

const double dt = 0.01;
 for( double t=0.0 ; t<10.0 ; t+= dt )
@@ -1179,7 +1175,6 @@
           with 4th order error estimation and coefficients introduced by Cash-Karp.
         

-

typedef runge_kutta_cash_karp54< state_type > error_stepper_type;
 
@@ -1192,7 +1187,6 @@ stepper exists which can be used via

-

typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
 controlled_stepper_type controlled_stepper;
@@ -1219,7 +1213,6 @@
           by using make_controlled:
         

-

integrate_adaptive( make_controlled< error_stepper_type >( 1.0e-10 , 1.0e-6 ) , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
 
@@ -1235,7 +1228,6 @@ is

-

integrate_adaptive( make_controlled( 1.0e-10 , 1.0e-6 , error_stepper_type() ) , harmonic_oscillator , x , 0.0 , 10.0 , 0.01 );
 
@@ -1251,7 +1243,6 @@ with all parameters is therefore

-

double abs_err = 1.0e-10 , rel_err = 1.0e-6 , a_x = 1.0 , a_dxdt = 1.0;
 controlled_stepper_type controlled_stepper( default_error_checker< double >( abs_err , rel_err , a_x , a_dxdt ) );
diff --git a/doc/boost_sandbox_numeric_odeint/tutorial/solar_system.html b/doc/boost_sandbox_numeric_odeint/tutorial/solar_system.html
index 05d75810..927dddd0 100644
--- a/doc/boost_sandbox_numeric_odeint/tutorial/solar_system.html
+++ b/doc/boost_sandbox_numeric_odeint/tutorial/solar_system.html
@@ -98,7 +98,6 @@
           space as well as the velocity. Therefore, we use the operators from Boost.Operators:
         

-

/*the point type */
 template< class T , size_t Dim >
@@ -136,7 +135,6 @@
           type we use std::tr1::array
         

-

// we simulate 5 planets and the sun
 const size_t n = 6;
@@ -159,7 +157,6 @@
           As system function we have to provide f(p) and f(q):
         

-

const double gravitational_constant = 2.95912208286e-4;
 
@@ -179,7 +176,6 @@
 

-

struct solar_system_momentum
 {
@@ -221,7 +217,6 @@
           apply here:
         

-

typedef symplectic_rkn_sb3a_mclachlan< container_type > stepper_type;
 const double dt = 100.0;
@@ -246,7 +241,6 @@
           The reference wrapper is also passed, but this is not a problem at all:
         

-

struct streaming_observer
 {
diff --git a/doc/boost_sandbox_numeric_odeint/tutorial/special_topics.html b/doc/boost_sandbox_numeric_odeint/tutorial/special_topics.html
index 774b2d1a..310cecac 100644
--- a/doc/boost_sandbox_numeric_odeint/tutorial/special_topics.html
+++ b/doc/boost_sandbox_numeric_odeint/tutorial/special_topics.html
@@ -63,7 +63,6 @@
           code is very simple
         

-

typedef boost::array< complex< double > , 1 > state_type;
 
@@ -89,7 +88,6 @@
           function. In this cast the Stuart-Landau oscillator looks like
         

-

double eta = 1.0;
 double alpha = 1.0;
@@ -111,7 +109,6 @@
           Integration is also very easy:
         

-

state_type x = {{ complex< double >( 1.0 , 0.0 ) }};
 
@@ -156,7 +153,6 @@
           components. Here, is how this function looks like
         

-

typedef vector< double > container_type;
 
@@ -208,7 +204,6 @@
           All this can be easily done with the following piece of code:
         

-

const size_t n = 64;
 container_type q( n , 0.0 ) , p( n , 0.0 );
@@ -235,7 +230,6 @@
           The observer is more or less trivial
         

-

struct streaming_observer
 {
@@ -304,7 +298,6 @@
           also need to store the individual frequencies of each oscillator.
         

-

typedef vector< double > container_type;
 
@@ -367,7 +360,6 @@
           and we record Z for different values of ε.
         

-

struct statistics_observer
 {
@@ -400,7 +392,6 @@
           2.
         

-

const size_t n = 16384;
 const double dt = 0.1;
@@ -458,7 +449,6 @@
           all quantities
         

-

// first increase fusion macro variables
 #define BOOST_FUSION_INVOKE_MAX_ARITY 15
@@ -500,7 +490,6 @@
           equivalent to the example in the first section of this tutorial
         

-

struct oscillator
 {
@@ -524,7 +513,6 @@
           Furthermore, the basic calculations are now performed by the fusion_algebra which must also be given.
         

-

typedef runge_kutta_dopri5< state_type , double , deriv_type , time_type , fusion_algebra > stepper_type;
 
@@ -539,7 +527,6 @@
           the observer is defined a bit different
         

-

struct streaming_observer
 {
@@ -611,7 +598,6 @@
           The definition of the system is
         

-

typedef boost::numeric::ublas::matrix< double > state_type;
 
@@ -705,7 +691,6 @@
           Here is a simple example:
         

-

typedef mpf_class value_type;
 typedef boost::array< value_type , 3 > state_type;
@@ -730,7 +715,6 @@
           which can be used
         

-

const int precision = 1024;
 mpf_set_default_prec( precision );
@@ -792,7 +776,6 @@
           is kept which should be changed whenever the states' size change.
         

-

typedef vector< double > coord_type;
 typedef pair< coord_type , coord_type > state_type;
@@ -878,7 +861,6 @@
           of 60.
         

-

//start with 60 sites
 const int N_start = 60;
@@ -914,7 +896,6 @@
           the start index of the potential changes in this case.
         

-

double t = 0.0;
 const double dt = 0.1;
@@ -947,7 +928,6 @@
           , p and distr.
         

-

void do_resize( coord_type &q , coord_type &p , coord_type &distr , const int N )
 {
diff --git a/doc/boost_sandbox_numeric_odeint/tutorial/stiff_systems.html b/doc/boost_sandbox_numeric_odeint/tutorial/stiff_systems.html
index 0d0fd902..a874b455 100644
--- a/doc/boost_sandbox_numeric_odeint/tutorial/stiff_systems.html
+++ b/doc/boost_sandbox_numeric_odeint/tutorial/stiff_systems.html
@@ -45,7 +45,6 @@
         is needed. Here is the definition of the above example
       

-

typedef boost::numeric::ublas::vector< double > vector_type;
 typedef boost::numeric::ublas::matrix< double > matrix_type;
@@ -82,7 +81,6 @@
         just templatize the operator():
       

-

typedef boost::numeric::ublas::vector< double > vector_type;
 typedef boost::numeric::ublas::matrix< double > matrix_type;
@@ -118,7 +116,6 @@
         all the other stepper:
       

-

vector_type x( 3 , 1.0 );
 
@@ -135,7 +132,6 @@
         5 method with step size control and dense output yields 1531 steps.
       

-

vector_type x2( 3 , 1.0 );
 
diff --git a/doc/boost_sandbox_numeric_odeint/tutorial/using_cuda_and_thrust.html b/doc/boost_sandbox_numeric_odeint/tutorial/using_cuda_and_thrust.html
index 7e40ed7d..ba716c69 100644
--- a/doc/boost_sandbox_numeric_odeint/tutorial/using_cuda_and_thrust.html
+++ b/doc/boost_sandbox_numeric_odeint/tutorial/using_cuda_and_thrust.html
@@ -87,7 +87,6 @@
           website.
         

-

//change this to float if your device does not support double computation
 typedef double value_type;
@@ -103,7 +102,7 @@
           a calculation on the GPU you usually have to call a global function like
           thrust::for_each, thrust::reduce,
           ... with an appropriate local functor which performs the basic operation.
-          An example is 
+          An example is
 

struct add_two
 {
@@ -138,7 +137,6 @@
           The mean field is calculated in a class mean_field_calculator
         

-

struct mean_field_calculator
 {
@@ -192,7 +190,6 @@
           GPU as well as on the CPU. The line
         

-

value_type sin_sum = thrust::reduce(
         thrust::make_transform_iterator( x.begin() , sin_functor() ) ,
@@ -208,7 +205,6 @@
           The system function is defined via
         

-

class phase_oscillator_ensemble
 {
@@ -264,7 +260,6 @@
           thrust algebra:
         

-

typedef runge_kutta4< state_type , value_type , state_type , value_type , thrust_algebra , thrust_operations > stepper_type;
 
@@ -274,7 +269,6 @@ Of course, you can also use a controlled or dense output stepper, e.g.

-

typedef runge_kutta_dopri5< state_type , value_type , state_type , value_type , thrust_algebra , thrust_operations > stepper_type;
 
@@ -285,14 +279,12 @@ an instance of the rhs class and using an integration function:

-

phase_oscillator_ensemble ensemble( N , 1.0 );
 

-

size_t steps1 = integrate_const( make_controlled( 1.0e-6 , 1.0e-6 , stepper_type() ) , boost::ref( ensemble ) , x , 0.0 , t_transients , dt );
 
@@ -334,7 +326,6 @@ in front. The full system class is:

-

//change this to host_vector< ... > of you want to run on CPU
 typedef thrust::device_vector< value_type > state_type;
@@ -423,7 +414,6 @@
           current state, hence the phase of each oscillator.
         

-

// create initial conditions and omegas on host:
 vector< value_type > x_host( N );
@@ -486,7 +476,6 @@
           >.
         

-

vector< value_type > beta_host( N );
 const value_type beta_min = 0.0 , beta_max = 56.0;
@@ -504,7 +493,6 @@
           calculates one particular realization of the Lorenz ensemble
         

-

struct lorenz_system
 {
@@ -586,7 +574,6 @@
           which is used within this observer is defined
         

-

struct lyap_functor
 {
@@ -626,7 +613,6 @@
           and written to stdout.
         

-

state_type x( 6 * N );
 
diff --git a/doc/index.html b/doc/index.html
index c83312c4..1a88e146 100644
--- a/doc/index.html
+++ b/doc/index.html
@@ -24,7 +24,7 @@
 
-

+

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)

@@ -125,7 +125,7 @@
- +

Last revised: March 29, 2012 at 12:51:03 GMT

Last revised: April 02, 2012 at 06:51:21 GMT


diff --git a/libs/numeric/odeint/doc/details_steppers.qbk b/libs/numeric/odeint/doc/details_steppers.qbk index 12bdff4e..8e37cd07 100644 --- a/libs/numeric/odeint/doc/details_steppers.qbk +++ b/libs/numeric/odeint/doc/details_steppers.qbk @@ -221,7 +221,7 @@ We will illustrate how to write your own stepper with the example of the stochas ['dx/dt = f(x) + g(x) __xi(t)] -where ['__xi] is Gaussian white noise with zero mean and a standard deviation ['__sigma]. ['f(x)] is said to be the deterministic part while [' g(x) __xi] 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 time step. BUt there exist many solvers for SDEs. A classical and easy method is the stochastic euler solver. It works by iterating +where ['__xi] is Gaussian white noise with zero mean and a standard deviation ['__sigma(t)]. ['f(x)] is said to be the deterministic part while [' g(x) __xi] 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 time step. BUt there exist many solvers for SDEs. A classical and easy method is the stochastic euler solver. It works by iterating ['x(t+__Delta t) = x(t) + __Delta t f(x(t)) + __Delta t[super 1/2] g(x) __xi(t)] @@ -231,17 +231,26 @@ Now we will implement this method. We will call the stepper `stochastic_euler`. [stochastic_euler_class_definition] -The types are needed in order to fullfil the stepper concept. An internal state is not needed in the stochastic euler, there are no temporaries needed during each step. So, in principle `state_type` and `deriv_type` are not needed here. Nevertheless, they must be present in order to fullfil the __stepper_concept and `std::vector< double >` is a good general candidate for the state and deriv type. +The types are needed in order to fulfill the stepper concept. As internal state and deriv type we use simple arrays in the stochastic euler, they are needed for the temporaries. The stepper has the order one which is returned from the `order()` function. +The system functions needs to calculate the deterministic and the stochastic part of our stochastic differential equation. So it might be suitable that the system function itself is a pair of functions, one for computing the 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 +[stochastic_euler_do_step] +This is all. It is quite simple and the stochastic euler stepper implement here is quite general. Of course it can be enhanced, for example +* use of operations and algebras as well as the resizing mechanism for maximal flexibility and portability +* use of `boost::ref` for the system functions +* use of `boost::range` for the state type in the `do_step` method +* ... -System 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 -steppers body +['dx/dt = - x + __xi] -do_step method +where __xi is Gaussian white noise with standard deviation ['__sigma]. + +[stochastic_euler_ornstein_uhlenbeck_det] examples, with integrate function diff --git a/libs/numeric/odeint/examples/stochastic_euler.cpp b/libs/numeric/odeint/examples/stochastic_euler.cpp index ca4987c3..8f995b5e 100644 --- a/libs/numeric/odeint/examples/stochastic_euler.cpp +++ b/libs/numeric/odeint/examples/stochastic_euler.cpp @@ -6,40 +6,59 @@ #include -struct sys -{ - template< class State > - void operator()( const State &x , State &dxdt ) - { - dxdt[0] = -x[0]; - } -}; - +/* //[ stochastic_euler_class_definition -class stochastic_euler +template< size_t N > class stochastic_euler { public: - typedef std::vector< double > state_type; - typedef std::vector< double > deriv_type; + typedef std::array< double , N > state_type; + typedef std::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 class stochastic_euler { public: - typedef std::vector< double > state_type; - typedef std::vector< double > deriv_type; + typedef std::array< double , N > state_type; + typedef std::array< double , N > deriv_type; typedef double value_type; typedef double time_type; typedef unsigned short order_type; @@ -48,32 +67,45 @@ public: static order_type order( void ) { return 1; } - template< class System , class State , class Time > - void do_step( System system , State &x , Time t , Time dt ) const + template< class System > + void do_step( System system , state_type &x , time_type t , time_type dt ) const { - State dxdt; - system.first( x , dxdt ); + deriv_type det , stoch ; + system.first( x , det ); + system.second( x , stoch ); for( size_t i=0 ; i -struct gen + + +//[ stochastic_euler_ornstein_uhlenbeck_def +const static size_t N = 1; +typedef std::array< double , N > state_type; + +struct ornstein_det { - Rng &rng; - Dist &dist; - gen( Rng &rng_ , Dist &dist_ ) : rng( rng_ ) , dist( dist_ ) { } - double operator()( void ) { return dist( rng ); } + void operator()( const state_type &x , state_type &dxdt ) const + { + dxdt[0] = -x[0]; + } }; -template< class Rng , class Dist > -gen< Rng , Dist > make_gen( Rng &rng , Dist &dist ) +struct ornstein_stoch { - return gen< Rng , Dist >( rng , dist ); -} + 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 ); + } +}; +//] struct streaming_observer { @@ -84,23 +116,14 @@ struct streaming_observer } }; -using namespace std; -using namespace boost::numeric::odeint; - -typedef std::array< double , 1 > state_type; - int main( int argc , char **argv ) { - - mt19937 rng; - normal_distribution<> dist( 0.0 , 1.0 ); - + using namespace std; + using namespace boost::numeric::odeint; double dt = 0.1; state_type x = {{ 1.0 }}; - streaming_observer obs; - cout << boost::is_same< streaming_observer , boost::numeric::odeint::null_observer >::value << endl; - integrate_const( stochastic_euler() , make_pair( sys() , make_gen( rng , dist ) ) , x , 0.0 , 10.0 , dt , obs ); + integrate_const( stochastic_euler< N >() , make_pair( ornstein_det() , ornstein_stoch( 1.0 ) ) , x , 0.0 , 10.0 , dt , streaming_observer() ); return 0; }