Squashed 'third_party/boostorg/odeint/' content from commit 6ff2719
Change-Id: If4892e29c1a5e6cf3a7aa51486a2725c251b0c7d
git-subtree-dir: third_party/boostorg/odeint
git-subtree-split: 6ff2719b6907b86596c3d43e88c1bcfdf29df560
diff --git a/examples/phase_oscillator_ensemble.cpp b/examples/phase_oscillator_ensemble.cpp
new file mode 100644
index 0000000..090fd58
--- /dev/null
+++ b/examples/phase_oscillator_ensemble.cpp
@@ -0,0 +1,151 @@
+/*
+ * 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;
+}