Merge commit '4ebe6417a5fce5f0994fc0c31bebf732be96a07c' as 'boost'

boost/libs/numeric/odeint/examples/chaotic_system.cpp

@@ -0,0 +1,119 @@
+/*
+ * 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. Furthermore,
+ * it shows how odeint interacts with boost.range.
+ *
+ * Copyright 2011-2012 Karsten Ahnert
+ * Copyright 2011-2013 Mario Mulansky
+ *
+ * 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)
+ */
+
+
+#include <iostream>
+#include <boost/array.hpp>
+
+#include <boost/numeric/odeint.hpp>
+
+#include "gram_schmidt.hpp"
+
+using namespace std;
+using namespace boost::numeric::odeint;
+
+
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+//[ system_function_without_perturbations
+struct lorenz
+{
+    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_ );
+
+        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 boost::array< double , N > state_type;
+typedef boost::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];
+    }
+}
+//]
+
+
+
+
+
+int main( int argc , char **argv )
+{
+    state_type x;
+    lyap_type lyap;
+
+    fill( x.begin() , x.end() , 0.0 );
+    x[0] = 10.0 ; x[1] = 10.0 ; x[2] = 5.0;
+
+    const double dt = 0.01;
+
+    //[ integrate_transients_with_range
+    // explicitly choose range_algebra to override default choice of array_algebra
+    runge_kutta4< state_type , double , state_type , double , range_algebra > rk4;
+
+    // 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 );
+
+    double t = 0.0;
+    size_t count = 0;
+    while( true )
+    {
+
+        t = integrate_n_steps( rk4 , lorenz_with_lyap , x , t , dt , 100 );
+        gram_schmidt< num_of_lyap >( x , lyap , n );
+        ++count;
+
+        if( !(count % 100000) )
+        {
+            cout << t;
+            for( size_t i=0 ; i<num_of_lyap ; ++i ) cout << "\t" << lyap[i] / t ;
+            cout << endl;
+        }
+    }
+    //]
+
+    return 0;
+}