In this example we present application of odeint to investigation of the properties of chaotic deterministic systems. In mathematical terms chaotic refers to an exponential growth of perturbations
where To calculate the Lyapunov exponents numerically one usually solves the equations of motion for To demonstrate how one can use odeint to determine the Lyapunov exponents we choose the Lorenz system. It is one of the most studied dynamical systems in the nonlinear dynamics community. For the standard parameters it possesses a strange attractor with non-integer dimension. The Lyapunov exponents take values of approximately 0.9, 0 and -12. The implementation of the Lorenz system is const double sigma = 10.0; const double R = 28.0; const double b = 8.0 / 3.0; typedef boost::array< double , 3 > lorenz_state_type; void lorenz( const lorenz_state_type &x , lorenz_state_type &dxdt , double t ) { 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]; } We need also to integrate the set of the perturbations. This is done in parallel to the original system, hence within one system function. Of course, we want to use the above definition of the Lorenz system, hence the definition of the system function including the Lorenz system itself and the perturbation could look like: const size_t n = 3; const size_t num_of_lyap = 3; const size_t N = n + n*num_of_lyap; typedef std::tr1::array< double , N > state_type; typedef std::tr1::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]; } } The perturbations are stored linearly in the struct lorenz { template< class StateIn , class StateOut , class Value > void operator()( const StateIn &x , StateOut &dxdt , Value t ) { 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]; } }; void lorenz_with_lyap( const state_type &x , state_type &dxdt , double t ) { lorenz()( x , dxdt , t ); ... } This works fine and state_type x; // initialize x.. explicit_rk4< state_type > rk4; integrate_n_steps( rk4 , lorenz_with_lyap , x , 0.0 , 0.01 , 1000 ); This code snippet performs 1000 steps with constant step size 0.01. A real world use case for the calculation of the Lyapunov exponents of Lorenz system would always include some transient steps, just to ensure that the current state lies on the attractor, hence it would look like state_type x; // initialize x explicit_rk4< state_type > rk4; integrate_n_steps( rk4 , lorenz , x , 0.0 , 0.01 , 1000 ); The problem is now, that void lorenz( double* x , double *dxdt , double t, void* params ) { ... } int system_length = 3; rk4( x , system_length , t , dt , lorenz ); But odeint supports a similar and much more sophisticated concept: Boost.Range. To make the steppers and the system ready to work with Boost.Range the system has to by changed: 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]; } }; This is in principle all. Now, we only have to call // 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 );
Having integrated a sufficient number of transients steps we are now able to calculate the Lyapunov exponents: - Initialize the perturbations. They are stored linearly behind the state of the Lorenz system. The perturbations are initialized such that
*p ij = δ ij*, where*p ij*is the*j*-component of the*i*.-th perturbation and*δ ij*is the Kronecker symbol. - Integrate 100 steps of the full system with perturbations
- Orthonormalize the perturbation using Gram-Schmidt orthonormalization algorithm.
- Repeat step 2 and 3. Every 10000 steps write the current Lyapunov exponent.
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; } } |

pBook > Solver:: Integrator:: C++ >