| /* |
| * phase_oscillator_ensemble.cpp |
| * |
| * Demonstrates the phase transition from an unsynchronized to an synchronized state. |
| * |
| * Copyright 2011-2012 Karsten Ahnert |
| * Copyright 2011-2012 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 <utility> |
| |
| #include <boost/numeric/odeint.hpp> |
| |
| #ifndef M_PI //not there on windows |
| #define M_PI 3.141592653589793 //... |
| #endif |
| |
| #include <boost/random.hpp> |
| |
| using namespace std; |
| using namespace boost::numeric::odeint; |
| |
| //[ phase_oscillator_ensemble_system_function |
| typedef vector< double > container_type; |
| |
| |
| pair< double , double > calc_mean_field( const container_type &x ) |
| { |
| size_t n = x.size(); |
| double cos_sum = 0.0 , sin_sum = 0.0; |
| for( size_t i=0 ; i<n ; ++i ) |
| { |
| cos_sum += cos( x[i] ); |
| sin_sum += sin( x[i] ); |
| } |
| cos_sum /= double( n ); |
| sin_sum /= double( n ); |
| |
| double K = sqrt( cos_sum * cos_sum + sin_sum * sin_sum ); |
| double Theta = atan2( sin_sum , cos_sum ); |
| |
| return make_pair( K , Theta ); |
| } |
| |
| |
| struct phase_ensemble |
| { |
| container_type m_omega; |
| double m_epsilon; |
| |
| phase_ensemble( const size_t n , double g = 1.0 , double epsilon = 1.0 ) |
| : m_omega( n , 0.0 ) , m_epsilon( epsilon ) |
| { |
| create_frequencies( g ); |
| } |
| |
| void create_frequencies( double g ) |
| { |
| boost::mt19937 rng; |
| boost::cauchy_distribution<> cauchy( 0.0 , g ); |
| boost::variate_generator< boost::mt19937&, boost::cauchy_distribution<> > gen( rng , cauchy ); |
| generate( m_omega.begin() , m_omega.end() , gen ); |
| } |
| |
| void set_epsilon( double epsilon ) { m_epsilon = epsilon; } |
| |
| double get_epsilon( void ) const { return m_epsilon; } |
| |
| void operator()( const container_type &x , container_type &dxdt , double /* t */ ) const |
| { |
| pair< double , double > mean = calc_mean_field( x ); |
| for( size_t i=0 ; i<x.size() ; ++i ) |
| dxdt[i] = m_omega[i] + m_epsilon * mean.first * sin( mean.second - x[i] ); |
| } |
| }; |
| //] |
| |
| |
| |
| //[ phase_oscillator_ensemble_observer |
| struct statistics_observer |
| { |
| double m_K_mean; |
| size_t m_count; |
| |
| statistics_observer( void ) |
| : m_K_mean( 0.0 ) , m_count( 0 ) { } |
| |
| template< class State > |
| void operator()( const State &x , double t ) |
| { |
| pair< double , double > mean = calc_mean_field( x ); |
| m_K_mean += mean.first; |
| ++m_count; |
| } |
| |
| double get_K_mean( void ) const { return ( m_count != 0 ) ? m_K_mean / double( m_count ) : 0.0 ; } |
| |
| void reset( void ) { m_K_mean = 0.0; m_count = 0; } |
| }; |
| //] |
| |
| |
| |
| |
| |
| |
| |
| |
| int main( int argc , char **argv ) |
| { |
| //[ phase_oscillator_ensemble_integration |
| const size_t n = 16384; |
| const double dt = 0.1; |
| |
| container_type x( n ); |
| |
| boost::mt19937 rng; |
| boost::uniform_real<> unif( 0.0 , 2.0 * M_PI ); |
| boost::variate_generator< boost::mt19937&, boost::uniform_real<> > gen( rng , unif ); |
| |
| // gamma = 1, the phase transition occurs at epsilon = 2 |
| phase_ensemble ensemble( n , 1.0 ); |
| statistics_observer obs; |
| |
| for( double epsilon = 0.0 ; epsilon < 5.0 ; epsilon += 0.1 ) |
| { |
| ensemble.set_epsilon( epsilon ); |
| obs.reset(); |
| |
| // start with random initial conditions |
| generate( x.begin() , x.end() , gen ); |
| |
| // calculate some transients steps |
| integrate_const( runge_kutta4< container_type >() , boost::ref( ensemble ) , x , 0.0 , 10.0 , dt ); |
| |
| // integrate and compute the statistics |
| integrate_const( runge_kutta4< container_type >() , boost::ref( ensemble ) , x , 0.0 , 100.0 , dt , boost::ref( obs ) ); |
| cout << epsilon << "\t" << obs.get_K_mean() << endl; |
| } |
| |
| |
| //] |
| |
| return 0; |
| } |