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