1 /*
2 * phase_oscillator_ensemble.cpp
3 *
4 * Demonstrates the phase transition from an unsynchronized to an synchronized state.
5 *
6 * Copyright 2011-2012 Karsten Ahnert
7 * Copyright 2011-2012 Mario Mulansky
8 * Distributed under the Boost Software License, Version 1.0. (See
9 * accompanying file LICENSE_1_0.txt or copy at
10 * http://www.boost.org/LICENSE_1_0.txt)
11 *
12 */
13
14 #include <iostream>
15 #include <utility>
16
17 #include <boost/numeric/odeint.hpp>
18
19 #ifndef M_PI //not there on windows
20 #define M_PI 3.141592653589793 //...
21 #endif
22
23 #include <boost/random.hpp>
24
25 using namespace std;
26 using namespace boost::numeric::odeint;
27
28 //[ phase_oscillator_ensemble_system_function
29 typedef vector< double > container_type;
30
31
calc_mean_field(const container_type & x)32 pair< double , double > calc_mean_field( const container_type &x )
33 {
34 size_t n = x.size();
35 double cos_sum = 0.0 , sin_sum = 0.0;
36 for( size_t i=0 ; i<n ; ++i )
37 {
38 cos_sum += cos( x[i] );
39 sin_sum += sin( x[i] );
40 }
41 cos_sum /= double( n );
42 sin_sum /= double( n );
43
44 double K = sqrt( cos_sum * cos_sum + sin_sum * sin_sum );
45 double Theta = atan2( sin_sum , cos_sum );
46
47 return make_pair( K , Theta );
48 }
49
50
51 struct phase_ensemble
52 {
53 container_type m_omega;
54 double m_epsilon;
55
phase_ensemblephase_ensemble56 phase_ensemble( const size_t n , double g = 1.0 , double epsilon = 1.0 )
57 : m_omega( n , 0.0 ) , m_epsilon( epsilon )
58 {
59 create_frequencies( g );
60 }
61
create_frequenciesphase_ensemble62 void create_frequencies( double g )
63 {
64 boost::mt19937 rng;
65 boost::cauchy_distribution<> cauchy( 0.0 , g );
66 boost::variate_generator< boost::mt19937&, boost::cauchy_distribution<> > gen( rng , cauchy );
67 generate( m_omega.begin() , m_omega.end() , gen );
68 }
69
set_epsilonphase_ensemble70 void set_epsilon( double epsilon ) { m_epsilon = epsilon; }
71
get_epsilonphase_ensemble72 double get_epsilon( void ) const { return m_epsilon; }
73
operator ()phase_ensemble74 void operator()( const container_type &x , container_type &dxdt , double /* t */ ) const
75 {
76 pair< double , double > mean = calc_mean_field( x );
77 for( size_t i=0 ; i<x.size() ; ++i )
78 dxdt[i] = m_omega[i] + m_epsilon * mean.first * sin( mean.second - x[i] );
79 }
80 };
81 //]
82
83
84
85 //[ phase_oscillator_ensemble_observer
86 struct statistics_observer
87 {
88 double m_K_mean;
89 size_t m_count;
90
statistics_observerstatistics_observer91 statistics_observer( void )
92 : m_K_mean( 0.0 ) , m_count( 0 ) { }
93
94 template< class State >
operator ()statistics_observer95 void operator()( const State &x , double t )
96 {
97 pair< double , double > mean = calc_mean_field( x );
98 m_K_mean += mean.first;
99 ++m_count;
100 }
101
get_K_meanstatistics_observer102 double get_K_mean( void ) const { return ( m_count != 0 ) ? m_K_mean / double( m_count ) : 0.0 ; }
103
resetstatistics_observer104 void reset( void ) { m_K_mean = 0.0; m_count = 0; }
105 };
106 //]
107
108
109
110
111
112
113
114
main(int argc,char ** argv)115 int main( int argc , char **argv )
116 {
117 //[ phase_oscillator_ensemble_integration
118 const size_t n = 16384;
119 const double dt = 0.1;
120
121 container_type x( n );
122
123 boost::mt19937 rng;
124 boost::uniform_real<> unif( 0.0 , 2.0 * M_PI );
125 boost::variate_generator< boost::mt19937&, boost::uniform_real<> > gen( rng , unif );
126
127 // gamma = 1, the phase transition occurs at epsilon = 2
128 phase_ensemble ensemble( n , 1.0 );
129 statistics_observer obs;
130
131 for( double epsilon = 0.0 ; epsilon < 5.0 ; epsilon += 0.1 )
132 {
133 ensemble.set_epsilon( epsilon );
134 obs.reset();
135
136 // start with random initial conditions
137 generate( x.begin() , x.end() , gen );
138
139 // calculate some transients steps
140 integrate_const( runge_kutta4< container_type >() , boost::ref( ensemble ) , x , 0.0 , 10.0 , dt );
141
142 // integrate and compute the statistics
143 integrate_const( runge_kutta4< container_type >() , boost::ref( ensemble ) , x , 0.0 , 100.0 , dt , boost::ref( obs ) );
144 cout << epsilon << "\t" << obs.get_K_mean() << endl;
145 }
146
147
148 //]
149
150 return 0;
151 }
152