Created
May 22, 2014 15:41
-
-
Save headmyshoulder/8a97787d63da8b8ffa32 to your computer and use it in GitHub Desktop.
lorenz system
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/* | |
* 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; | |
double R = 23.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 ) | |
{ | |
if( argc > 1 ) | |
{ | |
R = atof( argv[1] ); | |
} | |
cout << "R = " << R << endl; | |
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 << "\t" << x[0] << "\t" << x[1] << "\t" << x[2]; | |
cout << endl; | |
} | |
} | |
//] | |
return 0; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment