Boost-Commit :

Date view	Thread view	Subject view	Author view

Subject: [Boost-commit] svn:boost r71199 - in sandbox/odeint/branches/karsten/libs/numeric/odeint: doc doc/html doc/html/boost_sandbox_numeric_odeint examples
From: karsten.ahnert_at_[hidden]
Date: 2011-04-12 09:42:01

Author: karsten
Date: 2011-04-12 09:41:59 EDT (Tue, 12 Apr 2011)
New Revision: 71199
URL: http://svn.boost.org/trac/boost/changeset/71199

Log:
chaotic system example
Text files modified:
   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 | 114 ++++++++++++++++++++++++++++++---------
   sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/html/index.html | 2
   sandbox/odeint/branches/karsten/libs/numeric/odeint/doc/tutorial_chaotic_system.qbk | 34 +++++++++--
   sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/chaotic_system.cpp | 9 ++
   5 files changed, 122 insertions(+), 39 deletions(-)

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-12 09:41:59 EDT (Tue, 12 Apr 2011)
@@ -33,7 +33,7 @@
classes</a>
</h3></div></div></div>
<div class="table">
-<a name="id383118"></a>Table 1.4. Stepper Algorithms
+<a name="id383946"></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-12 09:41:59 EDT (Tue, 12 Apr 2011)
@@ -865,7 +865,7 @@
 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 [].
+ see for example .
 

 To calculate the Lyapunov exponents numerically one usually solves the equations
@@ -984,52 +984,110 @@

 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.
+ does not know that it should only use 3 element. In detail, odeint and its
+ steppers determine the length of the system under consideration be determing
+ the length of the state. In the classical solvers from the problem was solved
+ by pointer to the state and an appropriate length, something similar to
 

- hier gehts weiter!!!
+
+
+<pre class="programlisting">void lorenz( double* x , double *dxdt , double t, void* params )
+{
+ ...
+}
+
+int system_length = 3;
+rk4( x , system_length , t , dt , lorenz );
+</pre>
+
+ 
+
+ But odeint supports a similar and much more sophisticated concept: Boost.Range. To make the steppers and
+ the system ready to work with Boost.Range the system has to by changed:
 

 



-<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 )
+<pre class="programlisting">struct lorenz
{
- lorenz()( x , dxdt , t );
+ template< class State , class Deriv >
+ void operator()( const State &x_ , Deriv &dxdt_ , double t ) const
+ {
+ typename boost::range_iterator< const State >::type x = boost::begin( x_ );
+ typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ );

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

 

 

- show NR solution with pointers
+ This is in principle all. Now, we only have to call <code class="computeroutput">integrate_n_steps</code>
+ with an suited range:
 

- ranges
+ 
+
+
+
+<pre class="programlisting">// perform 10000 transient steps
+integrate_n_steps( rk4 , lorenz() , std::make_pair( x.begin() , x.begin() + n ) , 0.0 , dt , 10000 );
+</pre>
+
+ 
+
 

- templatized operator()
+ Having integrated a sufficient number of transients steps we are now able
+ to calculate the Lyapunov exponents:
+ 
+<div class="orderedlist"><ol type="1">
+<li>
+ First, we initialize the perturbations. They are stored linearly behind
+ the state of the Lorenz system. The perturbation are initialized such that
+ they fullfill i , pj > = delta ij where <,>
+ is the classical scalar product and delta ij is the Kronecker symbol.
+ </li>
+<li>
+ as
+ </li>
+</ol></div>
+
+ 
+
+
+
+<pre class="programlisting">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;
+fill( lyap.begin() , lyap.end() , 0.0 );
+
+double t = 0.0;
+size_t count = 0;
+while( true )
+{
+ t = integrate_n_steps( rk4 , lorenz_with_lyap , x , t , dt , 100 );
+ gram_schmitt( x , lyap , n , num_of_lyap );
+ ++count;
+
+ if( !(count % 100000) )
+ {
+ cout << t;
+ for( size_t i=0 ; i<num_of_lyap ; ++i ) cout << "\t" << lyap[i] / t ;
+ cout << endl;
+ }
+}
+</pre>
+
+ 
+
 

 gram schmitt
@@ -1362,7 +1420,7 @@
 The following table gives an overview over all examples.
 
<div class="table">
-<a name="id379970"></a>Table 1.3. Examples Overview
+<a name="id380798"></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-12 09:41:59 EDT (Tue, 12 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 06, 2011 at 19:34:18 GMT</td>
+<td align="left">Last revised: April 12, 2011 at 13:09:10 GMT</td>
<td align="right"><div class="copyright-footer"></div></td>
</tr></table>
<hr>

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-12 09:41:59 EDT (Tue, 12 Apr 2011)
@@ -1,10 +1,12 @@
[section Chaotic systems and Lyapunov exponents]

+[import ../examples/chaotic_system.cpp]
+
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 [].
+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.

@@ -87,22 +89,38 @@
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.
+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 detail, odeint and its steppers determine the length of the system under consideration be determing the length of the state. In the classical solvers from the problem was solved by pointer to the state and an appropriate length, something similar to
+
+``
+void lorenz( double* x , double *dxdt , double t, void* params )
+{
+ ...
+}
+
+int system_length = 3;
+rk4( x , system_length , t , dt , lorenz );
+``

-hier gehts weiter!!!
+But odeint supports a similar and much more sophisticated concept: __boost_range. To make the steppers and the system ready to work with __boost_range the system has to by changed:

+[system_function_without_perturbations]

-[import ../examples/chaotic_system.cpp]
-[system_function_with_perturbations]
+This is in principle all. Now, we only have to call `integrate_n_steps` with an suited range:
+
+[integrate_transients_with_range]
+
+Having integrated a sufficient number of transients steps we are now able to calculate the Lyapunov exponents:

+# First, we initialize the perturbations. They are stored linearly behind the state of the Lorenz system. The perturbation are initialized such that they fullfill
+[' = delta [subl ij] ]
+where <,> is the classical scalar product and delta [subl ij] is the Kronecker symbol.
+# as

+[lyapunov_full_code]

-show NR solution with pointers

-ranges

-templatized operator()

gram schmitt

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-12 09:41:59 EDT (Tue, 12 Apr 2011)
@@ -24,6 +24,7 @@
const double R = 28.0;
const double b = 8.0 / 3.0;

+//[ system_function_without_perturbations
struct lorenz
{
 template< class State , class Deriv >
@@ -37,6 +38,8 @@
 dxdt[2] = -b * x[2] + x[0] * x[1];
 }
};
+//]
+

//[ system_function_with_perturbations
@@ -77,9 +80,12 @@

 const double dt = 0.01;

- // 10000 transients steps
+ //[ integrate_transients_with_range
+ // perform 10000 transient steps
 integrate_n_steps( rk4 , lorenz() , std::make_pair( x.begin() , x.begin() + n ) , 0.0 , dt , 10000 );
+ //]

+ //[ lyapunov_full_code
 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;
 fill( lyap.begin() , lyap.end() , 0.0 );
@@ -99,6 +105,7 @@
 cout << endl;
 }
 }
+ //]

 return 0;
}

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