[Boost-commit] svn:boost r71008 - sandbox/odeint/branches/karsten/libs/numeric/odeint/examples

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