Boost logo

Boost-Commit :

Subject: [Boost-commit] svn:boost r71008 - sandbox/odeint/branches/karsten/libs/numeric/odeint/examples
From: karsten.ahnert_at_[hidden]
Date: 2011-04-05 14:31:35


Author: karsten
Date: 2011-04-05 14:31:34 EDT (Tue, 05 Apr 2011)
New Revision: 71008
URL: http://svn.boost.org/trac/boost/changeset/71008

Log:
chaotic system example finalized
Text files modified:
   sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/chaotic_system.cpp | 34 ++++----------------
   sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/gram_schmitt.hpp | 64 +++++++++++++++++++++------------------
   2 files changed, 41 insertions(+), 57 deletions(-)

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-05 14:31:34 EDT (Tue, 05 Apr 2011)
@@ -2,14 +2,14 @@
  * chaotic_system.cpp
  *
  * This example demonstrates how one can use odeint to determine the Lyapunov
- * exponents of a chaotic system namely the well known Lorenz system.
+ * exponents of a chaotic system namely the well known Lorenz system. Furthermore,
+ * it shows how odeint interacts with boost.range.
  *
  * Created on: Apr 5, 2011
  * Author: karsten
  */
 
 #include <iostream>
-
 #include <tr1/array>
 
 #include <boost/numeric/odeint.hpp>
@@ -56,14 +56,6 @@
     }
 }
 
-struct streaming
-{
- template< class State , class Time >
- void operator()( const State& x , const Time &t ) const
- {
- cout << t << "\t" << x[0] << "\t" << x[1] << "\t" << x[2] << "\n";
- }
-};
 
 int main( int argc , char **argv )
 {
@@ -74,7 +66,7 @@
     fill( x.begin() , x.end() , 0.0 );
     x[0] = 10.0 ; x[1] = 10.0 ; x[2] = 5.0;
 
- const double dt = 0.1;
+ const double dt = 0.01;
 
     // 10000 transients steps
     integrate_n_steps( rk4 , lorenz() , std::make_pair( x.begin() , x.begin() + n ) , 0.0 , dt , 10000 );
@@ -85,30 +77,18 @@
 
     double t = 0.0;
     size_t count = 0;
- for( size_t i=0 ; i<x.size() ; ++i ) cout << "\t" << x[i];
- cout << endl;
     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 == 1 )
- {
- for( size_t i=0 ; i<x.size() ; ++i ) cout << "\t" << x[i];
- cout << endl;
- }
+ gram_schmitt( x , lyap , n , num_of_lyap );
+ ++count;
 
-
- if( !(count % 1) )
+ if( !(count % 100000) )
         {
             cout << t;
- for( size_t i=0 ; i<num_of_lyap ; ++i )
- cout << "\t" << lyap[i] / t ;
- for( size_t i=0 ; i<x.size() ; ++i ) cout << "\t" << x[i];
+ for( size_t i=0 ; i<num_of_lyap ; ++i ) cout << "\t" << lyap[i] / t ;
             cout << endl;
         }
-
- if( count == 2 ) break;
     }
 
     return 0;

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/gram_schmitt.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/gram_schmitt.hpp (original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/examples/gram_schmitt.hpp 2011-04-05 14:31:34 EDT (Tue, 05 Apr 2011)
@@ -8,7 +8,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)
-*/
+ */
 
 #ifndef BOOST_NUMERIC_ODEINT_GRAM_SCHMITT_HPP_INCLUDED
 #define BOOST_NUMERIC_ODEINT_GRAM_SCHMITT_HPP_INCLUDED
@@ -21,27 +21,25 @@
 namespace numeric {
 namespace odeint {
 
- template< class Iterator , class T >
- void normalize( Iterator first , Iterator last , T norm )
- {
+template< class Iterator , class T >
+void normalize( Iterator first , Iterator last , T norm )
+{
         while( first != last ) *first++ /= norm;
- }
+}
 
- template< class Iterator , class T >
- void substract_vector( Iterator first1 , Iterator last1 ,
- Iterator first2 , T val )
- {
- while( first1 != last1 )
- *first1++ -= val * *first2++;
- }
-
- template< class StateType , class LyapType >
- void gram_schmitt( StateType &x , LyapType &lyap ,
- size_t n , size_t num_of_lyap )
- {
+template< class Iterator , class T >
+void substract_vector( Iterator first1 , Iterator last1 ,
+ Iterator first2 , T val )
+{
+ while( first1 != last1 ) *first1++ -= val * ( *first2++ );
+}
+
+template< class StateType , class LyapType >
+void gram_schmitt( StateType &x , LyapType &lyap , size_t n , size_t num_of_lyap )
+{
         if( !num_of_lyap ) return;
         if( ptrdiff_t( ( num_of_lyap + 1 ) * n ) != std::distance( x.begin() , x.end() ) )
- throw std::domain_error( "renormalization() : size of state does not match the number of lyapunov exponents." );
+ throw std::domain_error( "renormalization() : size of state does not match the number of lyapunov exponents." );
 
         typedef typename StateType::value_type value_type;
         typedef typename StateType::iterator iterator;
@@ -49,32 +47,38 @@
         value_type norm[num_of_lyap] , tmp[num_of_lyap];
         iterator first = x.begin() + n;
         iterator beg1 = first , end1 = first + n ;
-
+
         std::fill( norm , norm+num_of_lyap , 0.0 );
 
         // normalize first vector
         norm[0] = sqrt( std::inner_product( beg1 , end1 , beg1 , 0.0 ) );
         normalize( beg1 , end1 , norm[0] );
 
- beg1 += n ;
+ beg1 += n;
         end1 += n;
 
         for( size_t j=1 ; j<num_of_lyap ; ++j , beg1+=n , end1+=n )
         {
- for( size_t k=0 ; k<j ; ++k )
- tmp[k] = std::inner_product( beg1 , end1 , first + k*n , 0.0 );
+ for( size_t k=0 ; k<j ; ++k )
+ {
+ tmp[k] = std::inner_product( beg1 , end1 , first + k*n , 0.0 );
+// clog << j << " " << k << " " << tmp[k] << "\n";
+ }
+
+
+
+ for( size_t k=0 ; k<j ; ++k )
+ substract_vector( beg1 , end1 , first + k*n , tmp[k] );
 
- for( size_t k=0 ; k<j ; ++k )
- substract_vector( beg1 , end1 , first + k*n , tmp[k] );
-
- // nromalize j-th vector
- norm[j] = sqrt( std::inner_product( beg1 , end1 , beg1 , 0.0 ) );
- normalize( beg1 , end1 , norm[j] );
+ // nromalize j-th vector
+ norm[j] = sqrt( std::inner_product( beg1 , end1 , beg1 , 0.0 ) );
+// clog << j << " " << norm[j] << "\n";
+ normalize( beg1 , end1 , norm[j] );
         }
 
         for( size_t j=0 ; j<num_of_lyap ; j++ )
- lyap[j] += log( norm[j] );
- }
+ lyap[j] += log( norm[j] );
+}
 
 
 } // namespace odeint


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