1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Chaotic systems and Lyapunov exponents</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="solar_system.html" title="Solar system"> 10<link rel="next" href="stiff_systems.html" title="Stiff systems"> 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="solar_system.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="stiff_systems.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.chaotic_systems_and_lyapunov_exponents"></a><a class="link" href="chaotic_systems_and_lyapunov_exponents.html" title="Chaotic systems and Lyapunov exponents">Chaotic 28 systems and Lyapunov exponents</a> 29</h3></div></div></div> 30<p> 31 In this example we present application of odeint to investigation of the 32 properties of chaotic deterministic systems. In mathematical terms chaotic 33 refers to an exponential growth of perturbations <span class="emphasis"><em>δ x</em></span>. 34 In order to observe this exponential growth one usually solves the equations 35 for the tangential dynamics which is again an ordinary differential equation. 36 These equations are linear but time dependent and can be obtained via 37 </p> 38<p> 39 <span class="emphasis"><em>d δ x / dt = J(x) δ x</em></span> 40 </p> 41<p> 42 where <span class="emphasis"><em>J</em></span> is the Jacobian of the system under consideration. 43 <span class="emphasis"><em>δ x</em></span> can also be interpreted as a perturbation of the original 44 system. In principle <span class="emphasis"><em>n</em></span> of these perturbations exist, 45 they form a hypercube and evolve in the time. The Lyapunov exponents are 46 then defined as logarithmic growth rates of the perturbations. If one Lyapunov 47 exponent is larger then zero the nearby trajectories diverge exponentially 48 hence they are chaotic. If the largest Lyapunov exponent is zero one is usually 49 faced with periodic motion. In the case of a largest Lyapunov exponent smaller 50 then zero convergence to a fixed point is expected. More information's about 51 Lyapunov exponents and nonlinear dynamical systems can be found in many textbooks, 52 see for example: E. Ott "Chaos is Dynamical Systems", Cambridge. 53 </p> 54<p> 55 To calculate the Lyapunov exponents numerically one usually solves the equations 56 of motion for <span class="emphasis"><em>n</em></span> perturbations and orthonormalizes them 57 every <span class="emphasis"><em>k</em></span> steps. The Lyapunov exponent is the average 58 of the logarithm of the stretching factor of each perturbation. 59 </p> 60<p> 61 To demonstrate how one can use odeint to determine the Lyapunov exponents 62 we choose the Lorenz system. It is one of the most studied dynamical systems 63 in the nonlinear dynamics community. For the standard parameters it possesses 64 a strange attractor with non-integer dimension. The Lyapunov exponents take 65 values of approximately 0.9, 0 and -12. 66 </p> 67<p> 68 The implementation of the Lorenz system is 69 </p> 70<p> 71</p> 72<pre class="programlisting"><span class="keyword">const</span> <span class="keyword">double</span> <span class="identifier">sigma</span> <span class="special">=</span> <span class="number">10.0</span><span class="special">;</span> 73<span class="keyword">const</span> <span class="keyword">double</span> <span class="identifier">R</span> <span class="special">=</span> <span class="number">28.0</span><span class="special">;</span> 74<span class="keyword">const</span> <span class="keyword">double</span> <span class="identifier">b</span> <span class="special">=</span> <span class="number">8.0</span> <span class="special">/</span> <span class="number">3.0</span><span class="special">;</span> 75 76<span class="keyword">typedef</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">array</span><span class="special"><</span> <span class="keyword">double</span> <span class="special">,</span> <span class="number">3</span> <span class="special">></span> <span class="identifier">lorenz_state_type</span><span class="special">;</span> 77 78<span class="keyword">void</span> <span class="identifier">lorenz</span><span class="special">(</span> <span class="keyword">const</span> <span class="identifier">lorenz_state_type</span> <span class="special">&</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">lorenz_state_type</span> <span class="special">&</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">t</span> <span class="special">)</span> 79<span class="special">{</span> 80 <span class="identifier">dxdt</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">sigma</span> <span class="special">*</span> <span class="special">(</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">);</span> 81 <span class="identifier">dxdt</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">R</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">2</span><span class="special">];</span> 82 <span class="identifier">dxdt</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special">=</span> <span class="special">-</span><span class="identifier">b</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special">+</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">];</span> 83<span class="special">}</span> 84</pre> 85<p> 86 We need also to integrate the set of the perturbations. This is done in parallel 87 to the original system, hence within one system function. Of course, we want 88 to use the above definition of the Lorenz system, hence the definition of 89 the system function including the Lorenz system itself and the perturbation 90 could look like: 91 </p> 92<p> 93</p> 94<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">3</span><span class="special">;</span> 95<span class="keyword">const</span> <span class="identifier">size_t</span> <span class="identifier">num_of_lyap</span> <span class="special">=</span> <span class="number">3</span><span class="special">;</span> 96<span class="keyword">const</span> <span class="identifier">size_t</span> <span class="identifier">N</span> <span class="special">=</span> <span class="identifier">n</span> <span class="special">+</span> <span class="identifier">n</span><span class="special">*</span><span class="identifier">num_of_lyap</span><span class="special">;</span> 97 98<span class="keyword">typedef</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">tr1</span><span class="special">::</span><span class="identifier">array</span><span class="special"><</span> <span class="keyword">double</span> <span class="special">,</span> <span class="identifier">N</span> <span class="special">></span> <span class="identifier">state_type</span><span class="special">;</span> 99<span class="keyword">typedef</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">tr1</span><span class="special">::</span><span class="identifier">array</span><span class="special"><</span> <span class="keyword">double</span> <span class="special">,</span> <span class="identifier">num_of_lyap</span> <span class="special">></span> <span class="identifier">lyap_type</span><span class="special">;</span> 100 101<span class="keyword">void</span> <span class="identifier">lorenz_with_lyap</span><span class="special">(</span> <span class="keyword">const</span> <span class="identifier">state_type</span> <span class="special">&</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">state_type</span> <span class="special">&</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">t</span> <span class="special">)</span> 102<span class="special">{</span> 103 <span class="identifier">lorenz</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">,</span> <span class="identifier">dxdt</span> <span class="special">,</span> <span class="identifier">t</span> <span class="special">);</span> 104 105 <span class="keyword">for</span><span class="special">(</span> <span class="identifier">size_t</span> <span class="identifier">l</span><span class="special">=</span><span class="number">0</span> <span class="special">;</span> <span class="identifier">l</span><span class="special"><</span><span class="identifier">num_of_lyap</span> <span class="special">;</span> <span class="special">++</span><span class="identifier">l</span> <span class="special">)</span> 106 <span class="special">{</span> 107 <span class="keyword">const</span> <span class="keyword">double</span> <span class="special">*</span><span class="identifier">pert</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="number">3</span> <span class="special">+</span> <span class="identifier">l</span> <span class="special">*</span> <span class="number">3</span><span class="special">;</span> 108 <span class="keyword">double</span> <span class="special">*</span><span class="identifier">dpert</span> <span class="special">=</span> <span class="identifier">dxdt</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">+</span> <span class="number">3</span> <span class="special">+</span> <span class="identifier">l</span> <span class="special">*</span> <span class="number">3</span><span class="special">;</span> 109 <span class="identifier">dpert</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">=</span> <span class="special">-</span> <span class="identifier">sigma</span> <span class="special">*</span> <span class="identifier">pert</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">+</span> <span class="number">10.0</span> <span class="special">*</span> <span class="identifier">pert</span><span class="special">[</span><span class="number">1</span><span class="special">];</span> 110 <span class="identifier">dpert</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">=</span> <span class="special">(</span> <span class="identifier">R</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special">)</span> <span class="special">*</span> <span class="identifier">pert</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">pert</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">pert</span><span class="special">[</span><span class="number">2</span><span class="special">];</span> 111 <span class="identifier">dpert</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">pert</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">+</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">pert</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">b</span> <span class="special">*</span> <span class="identifier">pert</span><span class="special">[</span><span class="number">2</span><span class="special">];</span> 112 <span class="special">}</span> 113<span class="special">}</span> 114</pre> 115<p> 116 </p> 117<p> 118 The perturbations are stored linearly in the <code class="computeroutput"><span class="identifier">state_type</span></code> 119 behind the state of the Lorenz system. The problem of lorenz() and lorenz_with_lyap() having different 120 state types may be solved putting the Lorenz system inside a functor with 121 templatized arguments: 122 </p> 123<p> 124</p> 125<pre class="programlisting"><span class="keyword">struct</span> <span class="identifier">lorenz</span> 126<span class="special">{</span> 127 <span class="keyword">template</span><span class="special"><</span> <span class="keyword">class</span> <span class="identifier">StateIn</span> <span class="special">,</span> <span class="keyword">class</span> <span class="identifier">StateOut</span> <span class="special">,</span> <span class="keyword">class</span> <span class="identifier">Value</span> <span class="special">></span> 128 <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()(</span> <span class="keyword">const</span> <span class="identifier">StateIn</span> <span class="special">&</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">StateOut</span> <span class="special">&</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="identifier">Value</span> <span class="identifier">t</span> <span class="special">)</span> 129 <span class="special">{</span> 130 <span class="identifier">dxdt</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">sigma</span> <span class="special">*</span> <span class="special">(</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">);</span> 131 <span class="identifier">dxdt</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">R</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">2</span><span class="special">];</span> 132 <span class="identifier">dxdt</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special">=</span> <span class="special">-</span><span class="identifier">b</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special">+</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">];</span> 133 <span class="special">}</span> 134<span class="special">};</span> 135 136<span class="keyword">void</span> <span class="identifier">lorenz_with_lyap</span><span class="special">(</span> <span class="keyword">const</span> <span class="identifier">state_type</span> <span class="special">&</span><span class="identifier">x</span> <span class="special">,</span> <span class="identifier">state_type</span> <span class="special">&</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">t</span> <span class="special">)</span> 137<span class="special">{</span> 138 <span class="identifier">lorenz</span><span class="special">()(</span> <span class="identifier">x</span> <span class="special">,</span> <span class="identifier">dxdt</span> <span class="special">,</span> <span class="identifier">t</span> <span class="special">);</span> 139 <span class="special">...</span> 140<span class="special">}</span> 141</pre> 142<p> 143 This works fine and <code class="computeroutput"><span class="identifier">lorenz_with_lyap</span></code> 144 can be used for example via 145</p> 146<pre class="programlisting"><span class="identifier">state_type</span> <span class="identifier">x</span><span class="special">;</span> 147<span class="comment">// initialize x..</span> 148 149<span class="identifier">explicit_rk4</span><span class="special"><</span> <span class="identifier">state_type</span> <span class="special">></span> <span class="identifier">rk4</span><span class="special">;</span> 150<span class="identifier">integrate_n_steps</span><span class="special">(</span> <span class="identifier">rk4</span> <span class="special">,</span> <span class="identifier">lorenz_with_lyap</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">0.01</span> <span class="special">,</span> <span class="number">1000</span> <span class="special">);</span> 151</pre> 152<p> 153 This code snippet performs 1000 steps with constant step size 0.01. 154 </p> 155<p> 156 A real world use case for the calculation of the Lyapunov exponents of Lorenz 157 system would always include some transient steps, just to ensure that the 158 current state lies on the attractor, hence it would look like 159 </p> 160<p> 161</p> 162<pre class="programlisting"><span class="identifier">state_type</span> <span class="identifier">x</span><span class="special">;</span> 163<span class="comment">// initialize x</span> 164<span class="identifier">explicit_rk4</span><span class="special"><</span> <span class="identifier">state_type</span> <span class="special">></span> <span class="identifier">rk4</span><span class="special">;</span> 165<span class="identifier">integrate_n_steps</span><span class="special">(</span> <span class="identifier">rk4</span> <span class="special">,</span> <span class="identifier">lorenz</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">0.01</span> <span class="special">,</span> <span class="number">1000</span> <span class="special">);</span> 166</pre> 167<p> 168 The problem is now, that <code class="computeroutput"><span class="identifier">x</span></code> 169 is the full state containing also the perturbations and <code class="computeroutput"><span class="identifier">integrate_n_steps</span></code> 170 does not know that it should only use 3 elements. In detail, odeint and its 171 steppers determine the length of the system under consideration by determining 172 the length of the state. In the classical solvers, e.g. from Numerical Recipes, 173 the problem was solved by pointer to the state and an appropriate length, 174 something similar to 175 </p> 176<p> 177</p> 178<pre class="programlisting"><span class="keyword">void</span> <span class="identifier">lorenz</span><span class="special">(</span> <span class="keyword">double</span><span class="special">*</span> <span class="identifier">x</span> <span class="special">,</span> <span class="keyword">double</span> <span class="special">*</span><span class="identifier">dxdt</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">t</span><span class="special">,</span> <span class="keyword">void</span><span class="special">*</span> <span class="identifier">params</span> <span class="special">)</span> 179<span class="special">{</span> 180 <span class="special">...</span> 181<span class="special">}</span> 182 183<span class="keyword">int</span> <span class="identifier">system_length</span> <span class="special">=</span> <span class="number">3</span><span class="special">;</span> 184<span class="identifier">rk4</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">,</span> <span class="identifier">system_length</span> <span class="special">,</span> <span class="identifier">t</span> <span class="special">,</span> <span class="identifier">dt</span> <span class="special">,</span> <span class="identifier">lorenz</span> <span class="special">);</span> 185</pre> 186<p> 187 </p> 188<p> 189 But odeint supports a similar and much more sophisticated concept: <a href="http://www.boost.org/doc/libs/release/libs/range/" target="_top">Boost.Range</a>. 190 To make the steppers and the system ready to work with <a href="http://www.boost.org/doc/libs/release/libs/range/" target="_top">Boost.Range</a> 191 the system has to be changed: 192 </p> 193<p> 194</p> 195<pre class="programlisting"><span class="keyword">struct</span> <span class="identifier">lorenz</span> 196<span class="special">{</span> 197 <span class="keyword">template</span><span class="special"><</span> <span class="keyword">class</span> <span class="identifier">State</span> <span class="special">,</span> <span class="keyword">class</span> <span class="identifier">Deriv</span> <span class="special">></span> 198 <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="identifier">Deriv</span> <span class="special">&</span><span class="identifier">dxdt_</span> <span class="special">,</span> <span class="keyword">double</span> <span class="identifier">t</span> <span class="special">)</span> <span class="keyword">const</span> 199 <span class="special">{</span> 200 <span class="keyword">typename</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">range_iterator</span><span class="special"><</span> <span class="keyword">const</span> <span class="identifier">State</span> <span class="special">>::</span><span class="identifier">type</span> <span class="identifier">x</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">begin</span><span class="special">(</span> <span class="identifier">x_</span> <span class="special">);</span> 201 <span class="keyword">typename</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">range_iterator</span><span class="special"><</span> <span class="identifier">Deriv</span> <span class="special">>::</span><span class="identifier">type</span> <span class="identifier">dxdt</span> <span class="special">=</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">begin</span><span class="special">(</span> <span class="identifier">dxdt_</span> <span class="special">);</span> 202 203 <span class="identifier">dxdt</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">sigma</span> <span class="special">*</span> <span class="special">(</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">);</span> 204 <span class="identifier">dxdt</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">=</span> <span class="identifier">R</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">2</span><span class="special">];</span> 205 <span class="identifier">dxdt</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special">=</span> <span class="special">-</span><span class="identifier">b</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">2</span><span class="special">]</span> <span class="special">+</span> <span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">*</span> <span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">];</span> 206 <span class="special">}</span> 207<span class="special">};</span> 208</pre> 209<p> 210 </p> 211<p> 212 This is in principle all. Now, we only have to call <code class="computeroutput"><span class="identifier">integrate_n_steps</span></code> 213 with a range including only the first 3 components of <span class="emphasis"><em>x</em></span>: 214 </p> 215<p> 216</p> 217<pre class="programlisting"><span class="comment">// explicitly choose range_algebra to override default choice of array_algebra</span> 218<span class="identifier">runge_kutta4</span><span class="special"><</span> <span class="identifier">state_type</span> <span class="special">,</span> <span class="keyword">double</span> <span class="special">,</span> <span class="identifier">state_type</span> <span class="special">,</span> <span class="keyword">double</span> <span class="special">,</span> <span class="identifier">range_algebra</span> <span class="special">></span> <span class="identifier">rk4</span><span class="special">;</span> 219 220<span class="comment">// perform 10000 transient steps</span> 221<span class="identifier">integrate_n_steps</span><span class="special">(</span> <span class="identifier">rk4</span> <span class="special">,</span> <span class="identifier">lorenz</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">make_pair</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">begin</span><span class="special">()</span> <span class="special">+</span> <span class="identifier">n</span> <span class="special">)</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">,</span> <span class="identifier">dt</span> <span class="special">,</span> <span class="number">10000</span> <span class="special">);</span> 222</pre> 223<p> 224 </p> 225<div class="note"><table border="0" summary="Note"> 226<tr> 227<td rowspan="2" align="center" valign="top" width="25"><img alt="[Note]" src="../../../../../../../doc/src/images/note.png"></td> 228<th align="left">Note</th> 229</tr> 230<tr><td align="left" valign="top"><p> 231 Note that when using <a href="http://www.boost.org/doc/libs/release/libs/range/" target="_top">Boost.Range</a>, 232 we have to explicitly configure the stepper to use the <code class="computeroutput"><span class="identifier">range_algebra</span></code> 233 as otherwise odeint would automatically chose the <code class="computeroutput"><span class="identifier">array_algebra</span></code>, 234 which is incompatible with the usage of <a href="http://www.boost.org/doc/libs/release/libs/range/" target="_top">Boost.Range</a>, 235 because the original state_type is an <code class="computeroutput"><span class="identifier">array</span></code>. 236 </p></td></tr> 237</table></div> 238<p> 239 Having integrated a sufficient number of transients steps we are now able 240 to calculate the Lyapunov exponents: 241 </p> 242<div class="orderedlist"><ol class="orderedlist" type="1"> 243<li class="listitem"> 244 Initialize the perturbations. They are stored linearly behind the state 245 of the Lorenz system. The perturbations are initialized such that <span class="emphasis"><em>p 246 <sub>ij</sub> = δ <sub>ij</sub></em></span>, where <span class="emphasis"><em>p <sub>ij</sub></em></span> is the <span class="emphasis"><em>j</em></span>-component 247 of the <span class="emphasis"><em>i</em></span>.-th perturbation and <span class="emphasis"><em>δ <sub>ij</sub></em></span> 248 is the Kronecker symbol. 249 </li> 250<li class="listitem"> 251 Integrate 100 steps of the full system with perturbations 252 </li> 253<li class="listitem"> 254 Orthonormalize the perturbation using Gram-Schmidt orthonormalization 255 algorithm. 256 </li> 257<li class="listitem"> 258 Repeat step 2 and 3. Every 10000 steps write the current Lyapunov exponent. 259 </li> 260</ol></div> 261<p> 262</p> 263<pre class="programlisting"><span class="identifier">fill</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="identifier">n</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="number">0.0</span> <span class="special">);</span> 264<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">num_of_lyap</span> <span class="special">;</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">n</span><span class="special">+</span><span class="identifier">n</span><span class="special">*</span><span class="identifier">i</span><span class="special">+</span><span class="identifier">i</span><span class="special">]</span> <span class="special">=</span> <span class="number">1.0</span><span class="special">;</span> 265<span class="identifier">fill</span><span class="special">(</span> <span class="identifier">lyap</span><span class="special">.</span><span class="identifier">begin</span><span class="special">()</span> <span class="special">,</span> <span class="identifier">lyap</span><span class="special">.</span><span class="identifier">end</span><span class="special">()</span> <span class="special">,</span> <span class="number">0.0</span> <span class="special">);</span> 266 267<span class="keyword">double</span> <span class="identifier">t</span> <span class="special">=</span> <span class="number">0.0</span><span class="special">;</span> 268<span class="identifier">size_t</span> <span class="identifier">count</span> <span class="special">=</span> <span class="number">0</span><span class="special">;</span> 269<span class="keyword">while</span><span class="special">(</span> <span class="keyword">true</span> <span class="special">)</span> 270<span class="special">{</span> 271 272 <span class="identifier">t</span> <span class="special">=</span> <span class="identifier">integrate_n_steps</span><span class="special">(</span> <span class="identifier">rk4</span> <span class="special">,</span> <span class="identifier">lorenz_with_lyap</span> <span class="special">,</span> <span class="identifier">x</span> <span class="special">,</span> <span class="identifier">t</span> <span class="special">,</span> <span class="identifier">dt</span> <span class="special">,</span> <span class="number">100</span> <span class="special">);</span> 273 <span class="identifier">gram_schmidt</span><span class="special"><</span> <span class="identifier">num_of_lyap</span> <span class="special">>(</span> <span class="identifier">x</span> <span class="special">,</span> <span class="identifier">lyap</span> <span class="special">,</span> <span class="identifier">n</span> <span class="special">);</span> 274 <span class="special">++</span><span class="identifier">count</span><span class="special">;</span> 275 276 <span class="keyword">if</span><span class="special">(</span> <span class="special">!(</span><span class="identifier">count</span> <span class="special">%</span> <span class="number">100000</span><span class="special">)</span> <span class="special">)</span> 277 <span class="special">{</span> 278 <span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">t</span><span class="special">;</span> 279 <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">num_of_lyap</span> <span class="special">;</span> <span class="special">++</span><span class="identifier">i</span> <span class="special">)</span> <span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"\t"</span> <span class="special"><<</span> <span class="identifier">lyap</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">/</span> <span class="identifier">t</span> <span class="special">;</span> 280 <span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> 281 <span class="special">}</span> 282<span class="special">}</span> 283</pre> 284<p> 285 </p> 286<p> 287 The full code can be found here: <a href="https://github.com/headmyshoulder/odeint-v2/blob/master/examples/chaotic_system.cpp" target="_top">chaotic_system.cpp</a> 288 </p> 289</div> 290<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 291<td align="left"></td> 292<td align="right"><div class="copyright-footer">Copyright © 2009-2015 Karsten Ahnert and Mario Mulansky<p> 293 Distributed under the Boost Software License, Version 1.0. (See accompanying 294 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>) 295 </p> 296</div></td> 297</tr></table> 298<hr> 299<div class="spirit-nav"> 300<a accesskey="p" href="solar_system.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="stiff_systems.html"><img src="../../../../../../../doc/src/images/next.png" alt="Next"></a> 301</div> 302</body> 303</html> 304
