Еще одна важная высокомерная система соединенных обычных дифференциальных уравнений - это ансамбль N осцилляторов всех соединенных фаз [9] . Он определяется как

Природные частоты ω_​i каждого осциллятора следуют некоторому распределению и ε - сила соединения. Мы выбираем здесь дистрибутив Лоренца для ω_​i. Интересно, что фазовый переход можно наблюдать, если прочность соединения превышает критическое значение. Над этим значением задается синхронизация, и некоторые осцилляторы колебаются с той же частотой, несмотря на их различные естественные частоты. Переход также называется переходом Курамото. Его поведение можно проанализировать, используя среднее поле фазы

Определение функции системы сейчас немного сложнее, так как нам также нужно хранить отдельные частоты каждого осциллятора.

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

Z для различных значений ε. .

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

Теперь мы делаем несколько интеграций для различных значений ε и записи Z. Результат прекрасно подтверждает аналитический результат фазового перехода, т.е. в нашем примере стандартное отклонение Лоренца является 1 таким, что переход будет наблюдаться на ε = 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;
}