1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Ensembles of oscillators</title> 5<link rel="stylesheet" href="../../../../../../../doc/src/boostbook.css" type="text/css"> 6<meta name="generator" content="DocBook XSL Stylesheets V1.79.1"> 7<link rel="home" href="../../index.html" title="Chapter 1. Boost.Numeric.Odeint"> 8<link rel="up" href="../tutorial.html" title="Tutorial"> 9<link rel="prev" href="lattice_systems.html" title="Lattice systems"> 10<link rel="next" href="using_boost__units.html" title="Using boost::units"> 11</head> 12<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF"> 13<table cellpadding="2" width="100%"><tr> 14<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../logo.jpg"></td> 15<td align="center"><a href="../../../../../../../index.html">Home</a></td> 16<td align="center"><a href="../../../../../../../libs/libraries.htm">Libraries</a></td> 17<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td> 18<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td> 19<td align="center"><a href="../../../../../../../more/index.htm">More</a></td> 20</tr></table> 21<hr> 22<div class="spirit-nav"> 23<a accesskey="p" href="lattice_systems.html"><img src="../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../tutorial.html"><img src="../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="using_boost__units.html"><img src="../../../../../../../doc/src/images/next.png" alt="Next"></a> 24</div> 25<div class="section"> 26<div class="titlepage"><div><div><h3 class="title"> 27<a name="boost_numeric_odeint.tutorial.ensembles_of_oscillators"></a><a class="link" href="ensembles_of_oscillators.html" title="Ensembles of oscillators">Ensembles 28 of oscillators</a> 29</h3></div></div></div> 30<p> 31 Another important high dimensional system of coupled ordinary differential 32 equations is an ensemble of <span class="emphasis"><em>N</em></span> all-to-all coupled phase 33 oscillators <a class="link" href="../literature.html#synchronization_pikovsky_rosenblum">[9] </a>. 34 It is defined as 35 </p> 36<p> 37 <span class="emphasis"><em>dφ<sub>k</sub> / dt = ω<sub>k</sub> + ε / N Σ<sub>j</sub> sin( φ<sub>j</sub> - φ<sub>k</sub> )</em></span> 38 </p> 39<p> 40 The natural frequencies <span class="emphasis"><em>ω<sub>i</sub></em></span> of each oscillator follow 41 some distribution and <span class="emphasis"><em>ε</em></span> is the coupling strength. We 42 choose here a Lorentzian distribution for <span class="emphasis"><em>ω<sub>i</sub></em></span>. Interestingly 43 a phase transition can be observed if the coupling strength exceeds a critical 44 value. Above this value synchronization sets in and some of the oscillators 45 oscillate with the same frequency despite their different natural frequencies. 46 The transition is also called Kuramoto transition. Its behavior can be analyzed 47 by employing the mean field of the phase 48 </p> 49<p> 50 <span class="emphasis"><em>Z = K e<sup>i Θ</sup> = 1 / N Σ<sub>k</sub>e<sup>i φ<sub>k</sub></sup></em></span> 51 </p> 52<p> 53 The definition of the system function is now a bit more complex since we 54 also need to store the individual frequencies of each oscillator. 55 </p> 56<p> 57</p> 58<pre class="programlisting"><span class="keyword">typedef</span> <span class="identifier">vector</span><span class="special"><</span> <span class="keyword">double</span> <span class="special">></span> <span class="identifier">container_type</span><span class="special">;</span> 59 60 61<span class="identifier">pair</span><span class="special"><</span> <span class="keyword">double</span> <span class="special">,</span> <span class="keyword">double</span> <span class="special">></span> <span class="identifier">calc_mean_field</span><span class="special">(</span> <span class="keyword">const</span> <span class="identifier">container_type</span> <span class="special">&</span><span class="identifier">x</span> <span class="special">)</span> 62<span class="special">{</span> 63 <span class="identifier">size_t</span> <span class="identifier">n</span> <span class="special">=</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">size</span><span class="special">();</span> 64 <span class="keyword">double</span> <span class="identifier">cos_sum</span> <span class="special">=</span> <span class="number">0.0</span> <span class="special">,</span> <span class="identifier">sin_sum</span> <span class="special">=</span> <span class="number">0.0</span><span class="special">;</span> 65 <span class="keyword">for</span><span class="special">(</span> <span class="identifier">size_t</span> <span class="identifier">i</span><span class="special">=</span><span class="number">0</span> <span class="special">;</span> <span class="identifier">i</span><span class="special"><</span><span class="identifier">n</span> <span class="special">;</span> <span class="special">++</span><span class="identifier">i</span> <span class="special">)</span> 66 <span class="special">{</span> 67 <span class="identifier">cos_sum</span> <span class="special">+=</span> <span class="identifier">cos</span><span class="special">(</span> <span class="identifier">x</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">);</span> 68 <span class="identifier">sin_sum</span> <span class="special">+=</span> <span class="identifier">sin</span><span class="special">(</span> <span class="identifier">x</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">);</span> 69 <span class="special">}</span> 70 <span class="identifier">cos_sum</span> <span class="special">/=</span> <span class="keyword">double</span><span class="special">(</span> <span class="identifier">n</span> <span class="special">);</span> 71 <span class="identifier">sin_sum</span> <span class="special">/=</span> <span class="keyword">double</span><span class="special">(</span> <span class="identifier">n</span> <span class="special">);</span> 72 73 <span class="keyword">double</span> <span class="identifier">K</span> <span class="special">=</span> <span class="identifier">sqrt</span><span class="special">(</span> <span class="identifier">cos_sum</span> <span class="special">*</span> <span class="identifier">cos_sum</span> <span class="special">+</span> <span class="identifier">sin_sum</span> <span class="special">*</span> <span class="identifier">sin_sum</span> <span class="special">);</span> 74 <span class="keyword">double</span> <span class="identifier">Theta</span> <span class="special">=</span> <span class="identifier">atan2</span><span class="special">(</span> <span class="identifier">sin_sum</span> <span class="special">,</span> <span class="identifier">cos_sum</span> <span class="special">);</span> 75 76 <span class="keyword">return</span> <span class="identifier">make_pair</span><span class="special">(</span> <span class="identifier">K</span> <span class="special">,</span> <span class="identifier">Theta</span> <span class="special">);</span> 77<span class="special">}</span> 78 79 80<span class="keyword">struct</span> <span class="identifier">phase_ensemble</span> 81<span class="special">{</span> 82 <span class="identifier">container_type</span> <span class="identifier">m_omega</span><span class="special">;</span> 83 <span class="keyword">double</span> <span class="identifier">m_epsilon</span><span class="special">;</span> 84 85 <span class="identifier">phase_ensemble</span><span class="special">(</span> <span class="keyword">const</span> <span class="identifier">size_t</span> <span class="identifier">n</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">g</span> <span class="special">=</span> <span class="number">1.0</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">epsilon</span> <span class="special">=</span> <span class="number">1.0</span> <span class="special">)</span> 86 <span class="special">:</span> <span class="identifier">m_omega</span><span class="special">(</span> <span class="identifier">n</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">m_epsilon</span><span class="special">(</span> <span class="identifier">epsilon</span> <span class="special">)</span> 87 <span class="special">{</span> 88 <span class="identifier">create_frequencies</span><span class="special">(</span> <span class="identifier">g</span> <span class="special">);</span> 89 <span class="special">}</span> 90 91 <span class="keyword">void</span> <span class="identifier">create_frequencies</span><span class="special">(</span> <span class="keyword">double</span> <span class="identifier">g</span> <span class="special">)</span> 92 <span class="special">{</span> 93 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">mt19937</span> <span class="identifier">rng</span><span class="special">;</span> 94 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">cauchy_distribution</span><span class="special"><></span> <span class="identifier">cauchy</span><span class="special">(</span> <span class="number">0.0</span> <span class="special">,</span> <span class="identifier">g</span> <span class="special">);</span> 95 <span class="identifier">boost</span><span class="special">::</span><span class="identifier">variate_generator</span><span class="special"><</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">mt19937</span><span class="special">&,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">cauchy_distribution</span><span class="special"><></span> <span class="special">></span> <span class="identifier">gen</span><span class="special">(</span> <span class="identifier">rng</span> <span class="special">,</span> <span class="identifier">cauchy</span> <span class="special">);</span> 96 <span class="identifier">generate</span><span class="special">(</span> <span class="identifier">m_omega</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">m_omega</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">gen</span> <span class="special">);</span> 97 <span class="special">}</span> 98 99 <span class="keyword">void</span> <span class="identifier">set_epsilon</span><span class="special">(</span> <span class="keyword">double</span> <span class="identifier">epsilon</span> <span class="special">)</span> <span class="special">{</span> <span class="identifier">m_epsilon</span> <span class="special">=</span> <span class="identifier">epsilon</span><span class="special">;</span> <span class="special">}</span> 100 101 <span class="keyword">double</span> <span class="identifier">get_epsilon</span><span class="special">(</span> <span class="keyword">void</span> <span class="special">)</span> <span class="keyword">const</span> <span class="special">{</span> <span class="keyword">return</span> <span class="identifier">m_epsilon</span><span class="special">;</span> <span class="special">}</span> 102 103 <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="keyword">const</span> <span class="identifier">container_type</span> <span class="special">&</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">container_type</span> <span class="special">&</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="keyword">double</span> <span class="comment">/* t */</span> <span class="special">)</span> <span class="keyword">const</span> 104 <span class="special">{</span> 105 <span class="identifier">pair</span><span class="special"><</span> <span class="keyword">double</span> <span class="special">,</span> <span class="keyword">double</span> <span class="special">></span> <span class="identifier">mean</span> <span class="special">=</span> <span class="identifier">calc_mean_field</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">);</span> 106 <span class="keyword">for</span><span class="special">(</span> <span class="identifier">size_t</span> <span class="identifier">i</span><span class="special">=</span><span class="number">0</span> <span class="special">;</span> <span class="identifier">i</span><span class="special"><</span><span class="identifier">x</span><span class="special">.</span><span class="identifier">size</span><span class="special">()</span> <span class="special">;</span> <span class="special">++</span><span class="identifier">i</span> <span class="special">)</span> 107 <span class="identifier">dxdt</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">m_omega</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">+</span> <span class="identifier">m_epsilon</span> <span class="special">*</span> <span class="identifier">mean</span><span class="special">.</span><span class="identifier">first</span> <span class="special">*</span> <span class="identifier">sin</span><span class="special">(</span> <span class="identifier">mean</span><span class="special">.</span><span class="identifier">second</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">);</span> 108 <span class="special">}</span> 109<span class="special">};</span> 110</pre> 111<p> 112 </p> 113<p> 114 Note, that we have used <span class="emphasis"><em>Z</em></span> to simplify the equations 115 of motion. Next, we create an observer which computes the value of <span class="emphasis"><em>Z</em></span> 116 and we record <span class="emphasis"><em>Z</em></span> for different values of <span class="emphasis"><em>ε</em></span>. 117 </p> 118<p> 119</p> 120<pre class="programlisting"><span class="keyword">struct</span> <span class="identifier">statistics_observer</span> 121<span class="special">{</span> 122 <span class="keyword">double</span> <span class="identifier">m_K_mean</span><span class="special">;</span> 123 <span class="identifier">size_t</span> <span class="identifier">m_count</span><span class="special">;</span> 124 125 <span class="identifier">statistics_observer</span><span class="special">(</span> <span class="keyword">void</span> <span class="special">)</span> 126 <span class="special">:</span> <span class="identifier">m_K_mean</span><span class="special">(</span> <span class="number">0.0</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">m_count</span><span class="special">(</span> <span class="number">0</span> <span class="special">)</span> <span class="special">{</span> <span class="special">}</span> 127 128 <span class="keyword">template</span><span class="special"><</span> <span class="keyword">class</span> <span class="identifier">State</span> <span class="special">></span> 129 <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="keyword">const</span> <span class="identifier">State</span> <span class="special">&</span><span class="identifier">x</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">t</span> <span class="special">)</span> 130 <span class="special">{</span> 131 <span class="identifier">pair</span><span class="special"><</span> <span class="keyword">double</span> <span class="special">,</span> <span class="keyword">double</span> <span class="special">></span> <span class="identifier">mean</span> <span class="special">=</span> <span class="identifier">calc_mean_field</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">);</span> 132 <span class="identifier">m_K_mean</span> <span class="special">+=</span> <span class="identifier">mean</span><span class="special">.</span><span class="identifier">first</span><span class="special">;</span> 133 <span class="special">++</span><span class="identifier">m_count</span><span class="special">;</span> 134 <span class="special">}</span> 135 136 <span class="keyword">double</span> <span class="identifier">get_K_mean</span><span class="special">(</span> <span class="keyword">void</span> <span class="special">)</span> <span class="keyword">const</span> <span class="special">{</span> <span class="keyword">return</span> <span class="special">(</span> <span class="identifier">m_count</span> <span class="special">!=</span> <span class="number">0</span> <span class="special">)</span> <span class="special">?</span> <span class="identifier">m_K_mean</span> <span class="special">/</span> <span class="keyword">double</span><span class="special">(</span> <span class="identifier">m_count</span> <span class="special">)</span> <span class="special">:</span> <span class="number">0.0</span> <span class="special">;</span> <span class="special">}</span> 137 138 <span class="keyword">void</span> <span class="identifier">reset</span><span class="special">(</span> <span class="keyword">void</span> <span class="special">)</span> <span class="special">{</span> <span class="identifier">m_K_mean</span> <span class="special">=</span> <span class="number">0.0</span><span class="special">;</span> <span class="identifier">m_count</span> <span class="special">=</span> <span class="number">0</span><span class="special">;</span> <span class="special">}</span> 139<span class="special">};</span> 140</pre> 141<p> 142 </p> 143<p> 144 Now, we do several integrations for different values of <span class="emphasis"><em>ε</em></span> 145 and record <span class="emphasis"><em>Z</em></span>. The result nicely confirms the analytical 146 result of the phase transition, i.e. in our example the standard deviation 147 of the Lorentzian is 1 such that the transition will be observed at <span class="emphasis"><em>ε = 148 2</em></span>. 149 </p> 150<p> 151</p> 152<pre class="programlisting"><span class="keyword">const</span> <span class="identifier">size_t</span> <span class="identifier">n</span> <span class="special">=</span> <span class="number">16384</span><span class="special">;</span> 153<span class="keyword">const</span> <span class="keyword">double</span> <span class="identifier">dt</span> <span class="special">=</span> <span class="number">0.1</span><span class="special">;</span> 154 155<span class="identifier">container_type</span> <span class="identifier">x</span><span class="special">(</span> <span class="identifier">n</span> <span class="special">);</span> 156 157<span class="identifier">boost</span><span class="special">::</span><span class="identifier">mt19937</span> <span class="identifier">rng</span><span class="special">;</span> 158<span class="identifier">boost</span><span class="special">::</span><span class="identifier">uniform_real</span><span class="special"><></span> <span class="identifier">unif</span><span class="special">(</span> <span class="number">0.0</span> <span class="special">,</span> <span class="number">2.0</span> <span class="special">*</span> <span class="identifier">M_PI</span> <span class="special">);</span> 159<span class="identifier">boost</span><span class="special">::</span><span class="identifier">variate_generator</span><span class="special"><</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">mt19937</span><span class="special">&,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">uniform_real</span><span class="special"><></span> <span class="special">></span> <span class="identifier">gen</span><span class="special">(</span> <span class="identifier">rng</span> <span class="special">,</span> <span class="identifier">unif</span> <span class="special">);</span> 160 161<span class="comment">// gamma = 1, the phase transition occurs at epsilon = 2</span> 162<span class="identifier">phase_ensemble</span> <span class="identifier">ensemble</span><span class="special">(</span> <span class="identifier">n</span> <span class="special">,</span> <span class="number">1.0</span> <span class="special">);</span> 163<span class="identifier">statistics_observer</span> <span class="identifier">obs</span><span class="special">;</span> 164 165<span class="keyword">for</span><span class="special">(</span> <span class="keyword">double</span> <span class="identifier">epsilon</span> <span class="special">=</span> <span class="number">0.0</span> <span class="special">;</span> <span class="identifier">epsilon</span> <span class="special"><</span> <span class="number">5.0</span> <span class="special">;</span> <span class="identifier">epsilon</span> <span class="special">+=</span> <span class="number">0.1</span> <span class="special">)</span> 166<span class="special">{</span> 167 <span class="identifier">ensemble</span><span class="special">.</span><span class="identifier">set_epsilon</span><span class="special">(</span> <span class="identifier">epsilon</span> <span class="special">);</span> 168 <span class="identifier">obs</span><span class="special">.</span><span class="identifier">reset</span><span class="special">();</span> 169 170 <span class="comment">// start with random initial conditions</span> 171 <span class="identifier">generate</span><span class="special">(</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">x</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">gen</span> <span class="special">);</span> 172 173 <span class="comment">// calculate some transients steps</span> 174 <span class="identifier">integrate_const</span><span class="special">(</span> <span class="identifier">runge_kutta4</span><span class="special"><</span> <span class="identifier">container_type</span> <span class="special">>()</span> <span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">ref</span><span class="special">(</span> <span class="identifier">ensemble</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">x</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">,</span> <span class="number">10.0</span> <span class="special">,</span> <span class="identifier">dt</span> <span class="special">);</span> 175 176 <span class="comment">// integrate and compute the statistics</span> 177 <span class="identifier">integrate_const</span><span class="special">(</span> <span class="identifier">runge_kutta4</span><span class="special"><</span> <span class="identifier">container_type</span> <span class="special">>()</span> <span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">ref</span><span class="special">(</span> <span class="identifier">ensemble</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">x</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">,</span> <span class="number">100.0</span> <span class="special">,</span> <span class="identifier">dt</span> <span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">ref</span><span class="special">(</span> <span class="identifier">obs</span> <span class="special">)</span> <span class="special">);</span> 178 <span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">epsilon</span> <span class="special"><<</span> <span class="string">"\t"</span> <span class="special"><<</span> <span class="identifier">obs</span><span class="special">.</span><span class="identifier">get_K_mean</span><span class="special">()</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> 179<span class="special">}</span> 180</pre> 181<p> 182 </p> 183<p> 184 The full cpp file for this example can be found here <a href="https://github.com/headmyshoulder/odeint-v2/blob/master/examples/phase_oscillator_ensemble.cpp" target="_top">phase_oscillator_ensemble.cpp</a> 185 </p> 186</div> 187<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 188<td align="left"></td> 189<td align="right"><div class="copyright-footer">Copyright © 2009-2015 Karsten Ahnert and Mario Mulansky<p> 190 Distributed under the Boost Software License, Version 1.0. (See accompanying 191 file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>) 192 </p> 193</div></td> 194</tr></table> 195<hr> 196<div class="spirit-nav"> 197<a accesskey="p" href="lattice_systems.html"><img src="../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../tutorial.html"><img src="../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="using_boost__units.html"><img src="../../../../../../../doc/src/images/next.png" alt="Next"></a> 198</div> 199</body> 200</html> 201
