Boost-Commit :

Date view	Thread view	Subject view	Author view

Subject: [Boost-commit] svn:boost r71046 - in sandbox/odeint/branches/karsten: . libs/numeric/odeint/doc libs/numeric/odeint/doc/html libs/numeric/odeint/doc/html/boost_sandbox_numeric_odeint libs/numeric/odeint/examples
From: karsten.ahnert_at_[hidden]
Date: 2011-04-06 15:35:31

Next message: steven_at_[hidden]: "[Boost-commit] svn:boost r71047 - in trunk/libs/utility/doc: . html"
Previous message: erwann.rogard_at_[hidden]: "[Boost-commit] svn:boost r71045 - in sandbox/assign_v2/libs/assign/v2/test: option pipe/option"

Author: karsten
Date: 2011-04-06 15:35:28 EDT (Wed, 06 Apr 2011)
New Revision: 71046
URL: http://svn.boost.org/trac/boost/changeset/71046

Log:
chaotic system documentation
Text files modified:
   sandbox/odeint/branches/karsten/TODO | 1
   sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/html/boost_sandbox_numeric_odeint/reference.html | 2
   sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/html/boost_sandbox_numeric_odeint/tutorial.html | 183 +++++++++++++++++++++++++++++++++++++++
   sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/html/index.html | 2
   sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/odeint.qbk | 3
   sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial_chaotic_system.qbk | 98 +++++++++++++++++++++
   sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial_solar_system.qbk | 2
   sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/chaotic_system.cpp | 35 ++++--
   8 files changed, 307 insertions(+), 19 deletions(-)

Modified: sandbox/odeint/branches/karsten/TODO
==============================================================================
--- sandbox/odeint/branches/karsten/TODO (original)
+++ sandbox/odeint/branches/karsten/TODO 2011-04-06 15:35:28 EDT (Wed, 06 Apr 2011)
@@ -13,6 +13,7 @@
 * symplectic intergrators and ranges
* DOCUMENTATION
 * include links to the source of the examples
+ * use classical citations, via []
* GENERAL:
 * check if everywhere static_cast< value_type > is used
 * check header guards

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/html/boost_sandbox_numeric_odeint/reference.html
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/html/boost_sandbox_numeric_odeint/reference.html (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/html/boost_sandbox_numeric_odeint/reference.html 2011-04-06 15:35:28 EDT (Wed, 06 Apr 2011)
@@ -33,7 +33,7 @@
classes</a>
</h3></div></div></div>
<div class="table">
-<a name="id379450"></a>Table 1.4. Stepper Algorithms
+<a name="id383118"></a>Table 1.4. Stepper Algorithms
<div class="table-contents"><table class="table" summary="Stepper Algorithms">
<colgroup>
<col>

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/html/boost_sandbox_numeric_odeint/tutorial.html
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/html/boost_sandbox_numeric_odeint/tutorial.html (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/html/boost_sandbox_numeric_odeint/tutorial.html 2011-04-06 15:35:28 EDT (Wed, 06 Apr 2011)
@@ -600,7 +600,7 @@
 dqi / dt = pi
 

- dpi / dt = 1 / mi Σji Fij
+ dpboost_1_46_1.tari / dt = 1 / mi Σji Fij
 

 where pi is the momenta of object i.
@@ -844,6 +844,185 @@
 systems and Lyapunov exponents</a>
</h3></div></div></div>

+ Odeint can easily be used to investigate the properties of chaotic deterministic
+ systems. In mathematical terms chaotic refers to an exponential growth of
+ perturbations. In order to observe this exponential growth one usually solves
+ the equations for the tangential dynamics which is again an ordinary differential
+ equation. These equations are linear but time dependent and can be obtained
+ via
+ 
+
+ d delta x / dt = J(x) delta x
+ 
+
+ where J is the Jacobian of the system under consideration.
+ delta x cam also be interpreted as a perturbation of
+ the original system. In principle n of these perturbations
+ exit, they form a hypercube and evolve in the time. The Lyapunov exponents
+ are then defined as logarithmic growth rates of the perturbations. If one
+ Lyapunov exponent is larger then zero nearby trajectories divergy exponentially
+ hence they are chaotic. If the largest Lyapunov exponent is zero one is usually
+ faced with periodic motion. In the case of a largest Lyapunov exponent smaller
+ then zero the converges the a fixed point. More informations about Lyapunov
+ exponents and nonlinear dynamical systems can be found in many textbooks,
+ see for example [].
+ 
+
+ To calculate the Lyapunov exponents numerically one usually solves the equations
+ of motion for n perturbations und orthonormalized them
+ every k steps. The Lyapunov exponent is the average
+ the logarithm of the stretching factor of each perturbation.
+ 
+
+ To demonstrate how one can use odeint to determine the Lyapunov exponents
+ we choose the Lorenz system. It is one of the most studied dynamical systems
+ in the nonlinear dynamics community. For the standard parameters is possesses
+ a strange attractor with non-integer dimension. The Lyapunov exponents take
+ values of approximately 0.9, 0 and -12.
+ 
+
+ The implementation of the Lorenz system is
+ 
+
+
+
+<pre class="programlisting">const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+void lorenz( const state_type &x , state_type &dxdt , double t )
+{
+ dxdt[0] = sigma * ( x[1] - x[0] );
+ dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+ dxdt[2] = -b * x[2] + x[0] * x[1];
+}
+</pre>
+
+ Now, we need also need to integrate the set of the perturbations. This is
+ done in parallel to the original system, hence within one system function.
+ Of course, we want to use the above defintion of the Lorenz system, hence
+ the definition of the system function including the Lorenz system itself
+ and the perturbation could look like:
+ 
+
+
+
+<pre class="programlisting">const size_t n = 3;
+const size_t num_of_lyap = 3;
+const size_t N = n + n*num_of_lyap;
+
+typedef std::tr1::array< double , N > state_type;
+typedef std::tr1::array< double , num_of_lyap > lyap_type;
+
+void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t )
+{
+ lorenz( x , dxdt , t );
+
+ for( size_t l=0 ; l<num_of_lyap ; ++l )
+ {
+ const double *pert = x.begin() + 3 + l * 3;
+ double *dpert = dxdt.begin() + 3 + l * 3;
+ dpert[0] = - sigma * pert[0] + 10.0 * pert[1];
+ dpert[1] = ( R - x[2] ) * pert[0] - pert[1] - x[0] * pert[2];
+ dpert[2] = x[1] * pert[0] + x[0] * pert[1] - b * pert[2];
+ }
+}
+</pre>
+
+ 
+
+ But this approach will not work, since lorenz() and lorenz_with_lyap() have different state type. A simple
+ trick to over come this problem is put the Lorenz system inside a functor
+ with templatized arguments:
+ 
+
+
+
+<pre class="programlisting">struct lorenz
+{
+ template< class StateIn , class StateOut , class Value >
+ void operator()( const StateIn &x , StateOut &dxdt , Value t )
+ {
+ dxdt[0] = sigma * ( x[1] - x[0] );
+ dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+ dxdt[2] = -b * x[2] + x[0] * x[1];
+ }
+};
+
+void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t )
+{
+ lorenz()( x , dxdt , t );
+ ...
+}
+
+</pre>
+
+ This works fine and <code class="computeroutput">lorenz_with_lyap</code>
+ can be used for example via
+
+<pre class="programlisting">state_type x;
+// initialize x
+explicit_rk4< state_type > rk4;
+integrate_n_steps( rk4 , lorenz_with_lyap , x , 0.0 , 0.01 , 1000 );
+</pre>
+
+ This code snippet performs 1000 steps with constant step size 0.01.
+ 
+
+ A real world use case for the calculation of the Lyapunov exponents of Lorenz
+ system would always include some transient steps, just to ensure that the
+ current state lies on the attractor, hence it would look like
+ 
+
+
+
+<pre class="programlisting">state_type x;
+// initialize x
+explicit_rk4< state_type > rk4;
+integrate_n_steps( rk4 , lorenz , x , 0.0 , 0.01 , 1000 );
+</pre>
+
+ The problem is now, that <code class="computeroutput">x</code>
+ is the full state containing also the perturbations and <code class="computeroutput">integrate_n_steps</code>
+ does not know that it should only use 3 element. In classical solvers from
+ the Numerical recipies this problem was solved by pointers and an appropriate
+ length. But odeint supports a similar and much more sophisticated concept:
+ Boost.Range.
+ 
+
+ hier gehts weiter!!!
+ 
+
+ 
+
+
+
+<pre class="programlisting">const size_t n = 3;
+const size_t num_of_lyap = 3;
+const size_t N = n + n*num_of_lyap;
+
+typedef std::tr1::array< double , N > state_type;
+typedef std::tr1::array< double , num_of_lyap > lyap_type;
+
+void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t )
+{
+ lorenz()( x , dxdt , t );
+
+ for( size_t l=0 ; l<num_of_lyap ; ++l )
+ {
+ const double *pert = x.begin() + 3 + l * 3;
+ double *dpert = dxdt.begin() + 3 + l * 3;
+ dpert[0] = - sigma * pert[0] + 10.0 * pert[1];
+ dpert[1] = ( R - x[2] ) * pert[0] - pert[1] - x[0] * pert[2];
+ dpert[2] = x[1] * pert[0] + x[0] * pert[1] - b * pert[2];
+ }
+}
+</pre>
+
+ 
+
+ 
+
 show NR solution with pointers
 

@@ -1183,7 +1362,7 @@
 The following table gives an overview over all examples.
 
<div class="table">
-<a name="id376498"></a>Table 1.3. Examples Overview
+<a name="id379970"></a>Table 1.3. Examples Overview
<div class="table-contents"><table class="table" summary="Examples Overview">
<colgroup>
<col>

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/html/index.html
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/html/index.html (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/html/index.html 2011-04-06 15:35:28 EDT (Wed, 06 Apr 2011)
@@ -95,7 +95,7 @@
</div>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
-<td align="left">Last revised: April 05, 2011 at 07:39:00 GMT</td>
+<td align="left">Last revised: April 06, 2011 at 19:34:18 GMT</td>
<td align="right"><div class="copyright-footer"></div></td>
</tr></table>
<hr>

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/odeint.qbk
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/odeint.qbk (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/odeint.qbk 2011-04-06 15:35:28 EDT (Wed, 06 Apr 2011)
@@ -23,7 +23,8 @@
[@http://www.boost.org/doc/libs/release/doc/html/operators.html Boost.Operators]]
[def __boost_ref
[@http://www.boost.org/doc/libs/release/doc/html/ref.html `boost::ref`]]
-
+[def __boost_range
+ [@http://www.boost.org/doc/libs/release/libs/range/index.html `Boost.Range`]]

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial_chaotic_system.qbk
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial_chaotic_system.qbk (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial_chaotic_system.qbk 2011-04-06 15:35:28 EDT (Wed, 06 Apr 2011)
@@ -1,5 +1,103 @@
[section Chaotic systems and Lyapunov exponents]

+Odeint can easily be used to investigate the properties of chaotic deterministic systems. In mathematical terms chaotic refers to an exponential growth of perturbations. In order to observe this exponential growth one usually solves the equations for the tangential dynamics which is again an ordinary differential equation. These equations are linear but time dependent and can be obtained via
+
+['d delta x / dt = J(x) delta x]
+
+where ['J] is the Jacobian of the system under consideration. ['delta x] cam also be interpreted as a perturbation of the original system. In principle ['n] of these perturbations exit, they form a hypercube and evolve in the time. The Lyapunov exponents are then defined as logarithmic growth rates of the perturbations. If one Lyapunov exponent is larger then zero nearby trajectories divergy exponentially hence they are chaotic. If the largest Lyapunov exponent is zero one is usually faced with periodic motion. In the case of a largest Lyapunov exponent smaller then zero the converges the a fixed point. More informations about Lyapunov exponents and nonlinear dynamical systems can be found in many textbooks, see for example [].
+
+To calculate the Lyapunov exponents numerically one usually solves the equations of motion for ['n] perturbations und orthonormalized them every ['k] steps. The Lyapunov exponent is the average the logarithm of the stretching factor of each perturbation.
+
+To demonstrate how one can use odeint to determine the Lyapunov exponents we choose the Lorenz system. It is one of the most studied dynamical systems in the nonlinear dynamics community. For the standard parameters is possesses a strange attractor with non-integer dimension. The Lyapunov exponents take values of approximately 0.9, 0 and -12.
+
+The implementation of the Lorenz system is
+
+``
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+void lorenz( const state_type &x , state_type &dxdt , double t )
+{
+ dxdt[0] = sigma * ( x[1] - x[0] );
+ dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+ dxdt[2] = -b * x[2] + x[0] * x[1];
+}
+``
+Now, we need also need to integrate the set of the perturbations. This is done in parallel to the original system, hence within one system function. Of course, we want to use the above defintion of the Lorenz system, hence the definition of the system function including the Lorenz system itself and the perturbation could look like:
+
+``
+const size_t n = 3;
+const size_t num_of_lyap = 3;
+const size_t N = n + n*num_of_lyap;
+
+typedef std::tr1::array< double , N > state_type;
+typedef std::tr1::array< double , num_of_lyap > lyap_type;
+
+void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t )
+{
+ lorenz( x , dxdt , t );
+
+ for( size_t l=0 ; l<num_of_lyap ; ++l )
+ {
+ const double *pert = x.begin() + 3 + l * 3;
+ double *dpert = dxdt.begin() + 3 + l * 3;
+ dpert[0] = - sigma * pert[0] + 10.0 * pert[1];
+ dpert[1] = ( R - x[2] ) * pert[0] - pert[1] - x[0] * pert[2];
+ dpert[2] = x[1] * pert[0] + x[0] * pert[1] - b * pert[2];
+ }
+}
+``
+
+But this approach will not work, since '''lorenz()''' and '''lorenz_with_lyap()''' have different state type. A simple trick to over come this problem is put the Lorenz system inside a functor with templatized arguments:
+
+``
+struct lorenz
+{
+ template< class StateIn , class StateOut , class Value >
+ void operator()( const StateIn &x , StateOut &dxdt , Value t )
+ {
+ dxdt[0] = sigma * ( x[1] - x[0] );
+ dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+ dxdt[2] = -b * x[2] + x[0] * x[1];
+ }
+};
+
+void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t )
+{
+ lorenz()( x , dxdt , t );
+ ...
+}
+
+``
+This works fine and `lorenz_with_lyap` can be used for example via
+``
+state_type x;
+// initialize x
+explicit_rk4< state_type > rk4;
+integrate_n_steps( rk4 , lorenz_with_lyap , x , 0.0 , 0.01 , 1000 );
+``
+This code snippet performs 1000 steps with constant step size 0.01.
+
+A real world use case for the calculation of the Lyapunov exponents of Lorenz system would always include some transient steps, just to ensure that the current state lies on the attractor, hence it would look like
+
+``
+state_type x;
+// initialize x
+explicit_rk4< state_type > rk4;
+integrate_n_steps( rk4 , lorenz , x , 0.0 , 0.01 , 1000 );
+``
+The problem is now, that `x` is the full state containing also the perturbations and `integrate_n_steps` does not know that it should only use 3 element. In classical solvers from the Numerical recipies this problem was solved by pointers and an appropriate length. But odeint supports a similar and much more sophisticated concept: __boost_range.
+
+hier gehts weiter!!!
+
+
+[import ../examples/chaotic_system.cpp]
+[system_function_with_perturbations]
+
+
+
+
show NR solution with pointers

ranges

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial_solar_system.qbk
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial_solar_system.qbk (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial_solar_system.qbk 2011-04-06 15:35:28 EDT (Wed, 06 Apr 2011)
@@ -14,7 +14,7 @@

['dq[subl i] / dt = p[subl i]]

-['dp[subl i] / dt = 1 / m[subl i] __Sigma[subl ji] F[subl ij]]
+['dpboost_1_46_1.tar[subl i] / dt = 1 / m[subl i] __Sigma[subl ji] F[subl ij]]

where ['p[subl i]] is the momenta of object ['i]. The equations of motion can also be derived from the Hamiltonian

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/chaotic_system.cpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/chaotic_system.cpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/chaotic_system.cpp 2011-04-06 15:35:28 EDT (Wed, 06 Apr 2011)
@@ -19,14 +19,10 @@
using namespace std;
using namespace boost::numeric::odeint;

-const size_t n = 3;
-const size_t num_of_lyap = 3;
-const size_t N = n + n*num_of_lyap;
-
-typedef std::tr1::array< double , N > state_type;
-typedef std::tr1::array< double , num_of_lyap > lyap_type;
-const double dt = 0.01;

+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;

struct lorenz
{
@@ -36,12 +32,21 @@
 typename boost::range_iterator< const State >::type x = boost::begin( x_ );
 typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );

- dxdt[0]=10.0*(x[1]-x[0]);
- dxdt[1]=28.0*x[0]-x[1]-x[0]*x[2];
- dxdt[2]=-2.666666666666*x[2]+x[0]*x[1];
+ dxdt[0] = sigma * ( x[1] - x[0] );
+ dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+ dxdt[2] = -b * x[2] + x[0] * x[1];
 }
};

+
+//[ system_function_with_perturbations
+const size_t n = 3;
+const size_t num_of_lyap = 3;
+const size_t N = n + n*num_of_lyap;
+
+typedef std::tr1::array< double , N > state_type;
+typedef std::tr1::array< double , num_of_lyap > lyap_type;
+
void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t )
{
 lorenz()( x , dxdt , t );
@@ -50,11 +55,15 @@
 {
 const double *pert = x.begin() + 3 + l * 3;
 double *dpert = dxdt.begin() + 3 + l * 3;
- dpert[0] = -10.0 * pert[0] + 10.0 * pert[1];
- dpert[1] = ( 28.0 - x[2] ) * pert[0] - pert[1] - x[0] * pert[2];
- dpert[2] = x[1] * pert[0] + x[0] * pert[1] - 2.666666666666 * pert[2];
+ dpert[0] = - sigma * pert[0] + 10.0 * pert[1];
+ dpert[1] = ( R - x[2] ) * pert[0] - pert[1] - x[0] * pert[2];
+ dpert[2] = x[1] * pert[0] + x[0] * pert[1] - b * pert[2];
 }
}
+//]
+
+
+

int main( int argc , char **argv )

Next message: steven_at_[hidden]: "[Boost-commit] svn:boost r71047 - in trunk/libs/utility/doc: . html"
Previous message: erwann.rogard_at_[hidden]: "[Boost-commit] svn:boost r71045 - in sandbox/assign_v2/libs/assign/v2/test: option pipe/option"

Date view	Thread view	Subject view	Author view

Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk