1 /* 2 * chaotic_system.cpp 3 * 4 * This example demonstrates how one can use odeint to determine the Lyapunov 5 * exponents of a chaotic system namely the well known Lorenz system. Furthermore, 6 * it shows how odeint interacts with boost.range. 7 * 8 * Copyright 2011-2012 Karsten Ahnert 9 * Copyright 2011-2013 Mario Mulansky 10 * 11 * Distributed under the Boost Software License, Version 1.0. 12 * (See accompanying file LICENSE_1_0.txt or 13 * copy at http://www.boost.org/LICENSE_1_0.txt) 14 */ 15 16 17 #include <iostream> 18 #include <boost/array.hpp> 19 20 #include <boost/numeric/odeint.hpp> 21 22 #include "gram_schmidt.hpp" 23 24 using namespace std; 25 using namespace boost::numeric::odeint; 26 27 28 const double sigma = 10.0; 29 const double R = 28.0; 30 const double b = 8.0 / 3.0; 31 32 //[ system_function_without_perturbations 33 struct lorenz 34 { 35 template< class State , class Deriv > operator ()lorenz36 void operator()( const State &x_ , Deriv &dxdt_ , double t ) const 37 { 38 typename boost::range_iterator< const State >::type x = boost::begin( x_ ); 39 typename boost::range_iterator< Deriv >::type dxdt = boost::begin( dxdt_ ); 40 41 dxdt[0] = sigma * ( x[1] - x[0] ); 42 dxdt[1] = R * x[0] - x[1] - x[0] * x[2]; 43 dxdt[2] = -b * x[2] + x[0] * x[1]; 44 } 45 }; 46 //] 47 48 49 50 //[ system_function_with_perturbations 51 const size_t n = 3; 52 const size_t num_of_lyap = 3; 53 const size_t N = n + n*num_of_lyap; 54 55 typedef boost::array< double , N > state_type; 56 typedef boost::array< double , num_of_lyap > lyap_type; 57 lorenz_with_lyap(const state_type & x,state_type & dxdt,double t)58 void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t ) 59 { 60 lorenz()( x , dxdt , t ); 61 62 for( size_t l=0 ; l<num_of_lyap ; ++l ) 63 { 64 const double *pert = x.begin() + 3 + l * 3; 65 double *dpert = dxdt.begin() + 3 + l * 3; 66 dpert[0] = - sigma * pert[0] + 10.0 * pert[1]; 67 dpert[1] = ( R - x[2] ) * pert[0] - pert[1] - x[0] * pert[2]; 68 dpert[2] = x[1] * pert[0] + x[0] * pert[1] - b * pert[2]; 69 } 70 } 71 //] 72 73 74 75 76 main(int argc,char ** argv)77 int main( int argc , char **argv ) 78 { 79 state_type x; 80 lyap_type lyap; 81 82 fill( x.begin() , x.end() , 0.0 ); 83 x[0] = 10.0 ; x[1] = 10.0 ; x[2] = 5.0; 84 85 const double dt = 0.01; 86 87 //[ integrate_transients_with_range 88 // explicitly choose range_algebra to override default choice of array_algebra 89 runge_kutta4< state_type , double , state_type , double , range_algebra > rk4; 90 91 // perform 10000 transient steps 92 integrate_n_steps( rk4 , lorenz() , std::make_pair( x.begin() , x.begin() + n ) , 0.0 , dt , 10000 ); 93 //] 94 95 //[ lyapunov_full_code 96 fill( x.begin()+n , x.end() , 0.0 ); 97 for( size_t i=0 ; i<num_of_lyap ; ++i ) x[n+n*i+i] = 1.0; 98 fill( lyap.begin() , lyap.end() , 0.0 ); 99 100 double t = 0.0; 101 size_t count = 0; 102 while( true ) 103 { 104 105 t = integrate_n_steps( rk4 , lorenz_with_lyap , x , t , dt , 100 ); 106 gram_schmidt< num_of_lyap >( x , lyap , n ); 107 ++count; 108 109 if( !(count % 100000) ) 110 { 111 cout << t; 112 for( size_t i=0 ; i<num_of_lyap ; ++i ) cout << "\t" << lyap[i] / t ; 113 cout << endl; 114 } 115 } 116 //] 117 118 return 0; 119 } 120