Another important high dimensional system of coupled ordinary differential equations is an ensemble of N all-to-all coupled phase oscillators [9] . It is defined as

dφ_​k / dt = ω_​k + ε / N Σ_​j sin( φ_​j - φ_​k )

The natural frequencies ω_​i of each oscillator follow some distribution and ε is the coupling strength. We choose here a Lorentzian distribution for ω_​i. Interestingly a phase transition can be observed if the coupling strength exceeds a critical value. Above this value synchronization sets in and some of the oscillators oscillate with the same frequency despite their different natural frequencies. The transition is also called Kuramoto transition. Its behavior can be analyzed by employing the mean field of the phase

Z = K e^{i Θ} = 1 / N Σ_​ke^{i φ_​k}

The definition of the system function is now a bit more complex since we also need to store the individual frequencies of each oscillator.

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] );
    }
};

Note, that we have used Z to simplify the equations of motion. Next, we create an observer which computes the value of Z and we record Z for different values of ε.

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; }
};

Now, we do several integrations for different values of ε and record Z. The result nicely confirms the analytical result of the phase transition, i.e. in our example the standard deviation of the Lorentzian is 1 such that the transition will be observed at ε = 2.

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;
}

The full cpp file for this example can be found here phase_oscillator_ensemble.cpp