152 lines
3.7 KiB
C++
152 lines
3.7 KiB
C++
/*
|
|
* 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;
|
|
}
|