1<html> 2<head> 3<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> 4<title>Short Example</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="../getting_started.html" title="Getting started"> 9<link rel="prev" href="usage__compilation__headers.html" title="Usage, Compilation, Headers"> 10<link rel="next" href="../tutorial.html" title="Tutorial"> 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="usage__compilation__headers.html"><img src="../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../getting_started.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="../tutorial.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.getting_started.short_example"></a><a class="link" href="short_example.html" title="Short Example">Short 28 Example</a> 29</h3></div></div></div> 30<p> 31 Imaging, you want to numerically integrate a harmonic oscillator with friction. 32 The equations of motion are given by <span class="emphasis"><em>x'' = -x + γ x'</em></span>. 33 Odeint only deals with first order ODEs that have no higher derivatives than 34 x' involved. However, any higher order ODE can be transformed to a system 35 of first order ODEs by introducing the new variables <span class="emphasis"><em>q=x</em></span> 36 and <span class="emphasis"><em>p=x'</em></span> such that <span class="emphasis"><em>w=(q,p)</em></span>. To 37 apply numerical integration one first has to design the right hand side of 38 the equation <span class="emphasis"><em>w' = f(w) = (p,-q+γ p)</em></span>: 39 </p> 40<p> 41</p> 42<pre class="programlisting"><span class="comment">/* The type of container used to hold the state vector */</span> 43<span class="keyword">typedef</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span> <span class="keyword">double</span> <span class="special">></span> <span class="identifier">state_type</span><span class="special">;</span> 44 45<span class="keyword">const</span> <span class="keyword">double</span> <span class="identifier">gam</span> <span class="special">=</span> <span class="number">0.15</span><span class="special">;</span> 46 47<span class="comment">/* The rhs of x' = f(x) */</span> 48<span class="keyword">void</span> <span class="identifier">harmonic_oscillator</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">const</span> <span class="keyword">double</span> <span class="comment">/* t */</span> <span class="special">)</span> 49<span class="special">{</span> 50 <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">x</span><span class="special">[</span><span class="number">1</span><span class="special">];</span> 51 <span class="identifier">dxdt</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">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">gam</span><span class="special">*</span><span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">];</span> 52<span class="special">}</span> 53</pre> 54<p> 55 </p> 56<p> 57 Here we chose <code class="computeroutput"><span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span></code> 58 as the state type, but others are also possible, for example <code class="computeroutput"><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">2</span><span class="special">></span></code>. odeint is designed in such a way that 59 you can easily use your own state types. Next, the ODE is defined which is 60 in this case a simple function calculating <span class="emphasis"><em>f(x)</em></span>. The 61 parameter signature of this function is crucial: the integration methods 62 will always call them in the form <code class="computeroutput"><span class="identifier">f</span><span class="special">(</span><span class="identifier">x</span><span class="special">,</span> 63 <span class="identifier">dxdt</span><span class="special">,</span> 64 <span class="identifier">t</span><span class="special">)</span></code> 65 (there are exceptions for some special routines). So, even if there is no 66 explicit time dependence, one has to define <code class="computeroutput"><span class="identifier">t</span></code> 67 as a function parameter. 68 </p> 69<p> 70 Now, we have to define the initial state from which the integration should 71 start: 72 </p> 73<p> 74</p> 75<pre class="programlisting"><span class="identifier">state_type</span> <span class="identifier">x</span><span class="special">(</span><span class="number">2</span><span class="special">);</span> 76<span class="identifier">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">=</span> <span class="number">1.0</span><span class="special">;</span> <span class="comment">// start at x=1.0, p=0.0</span> 77<span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="special">=</span> <span class="number">0.0</span><span class="special">;</span> 78</pre> 79<p> 80 </p> 81<p> 82 For the integration itself we'll use the <code class="computeroutput"><span class="identifier">integrate</span></code> 83 function, which is a convenient way to get quick results. It is based on 84 the error-controlled <code class="computeroutput"><span class="identifier">runge_kutta54_cash_karp</span></code> 85 stepper (5th order) and uses adaptive step-size. 86 </p> 87<p> 88</p> 89<pre class="programlisting"><span class="identifier">size_t</span> <span class="identifier">steps</span> <span class="special">=</span> <span class="identifier">integrate</span><span class="special">(</span> <span class="identifier">harmonic_oscillator</span> <span class="special">,</span> 90 <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="number">0.1</span> <span class="special">);</span> 91</pre> 92<p> 93 </p> 94<p> 95 The integrate function expects as parameters the rhs of the ode as defined 96 above, the initial state <code class="computeroutput"><span class="identifier">x</span></code>, 97 the start-and end-time of the integration as well as the initial time step=size. 98 Note, that <code class="computeroutput"><span class="identifier">integrate</span></code> uses 99 an adaptive step-size during the integration steps so the time points will 100 not be equally spaced. The integration returns the number of steps that were 101 applied and updates x which is set to the approximate solution of the ODE 102 at the end of integration. 103 </p> 104<p> 105 It is also possible to represent the ode system as a class. The rhs must 106 then be implemented as a functor - a class with an overloaded function call 107 operator: 108 </p> 109<p> 110</p> 111<pre class="programlisting"><span class="comment">/* The rhs of x' = f(x) defined as a class */</span> 112<span class="keyword">class</span> <span class="identifier">harm_osc</span> <span class="special">{</span> 113 114 <span class="keyword">double</span> <span class="identifier">m_gam</span><span class="special">;</span> 115 116<span class="keyword">public</span><span class="special">:</span> 117 <span class="identifier">harm_osc</span><span class="special">(</span> <span class="keyword">double</span> <span class="identifier">gam</span> <span class="special">)</span> <span class="special">:</span> <span class="identifier">m_gam</span><span class="special">(</span><span class="identifier">gam</span><span class="special">)</span> <span class="special">{</span> <span class="special">}</span> 118 119 <span class="keyword">void</span> <span class="keyword">operator</span><span class="special">()</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">const</span> <span class="keyword">double</span> <span class="comment">/* t */</span> <span class="special">)</span> 120 <span class="special">{</span> 121 <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">x</span><span class="special">[</span><span class="number">1</span><span class="special">];</span> 122 <span class="identifier">dxdt</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">x</span><span class="special">[</span><span class="number">0</span><span class="special">]</span> <span class="special">-</span> <span class="identifier">m_gam</span><span class="special">*</span><span class="identifier">x</span><span class="special">[</span><span class="number">1</span><span class="special">];</span> 123 <span class="special">}</span> 124<span class="special">};</span> 125</pre> 126<p> 127 </p> 128<p> 129 which can be used via 130 </p> 131<p> 132</p> 133<pre class="programlisting"><span class="identifier">harm_osc</span> <span class="identifier">ho</span><span class="special">(</span><span class="number">0.15</span><span class="special">);</span> 134<span class="identifier">steps</span> <span class="special">=</span> <span class="identifier">integrate</span><span class="special">(</span> <span class="identifier">ho</span> <span class="special">,</span> 135 <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="number">0.1</span> <span class="special">);</span> 136</pre> 137<p> 138 </p> 139<p> 140 In order to observe the solution during the integration steps all you have 141 to do is to provide a reasonable observer. An example is 142 </p> 143<p> 144</p> 145<pre class="programlisting"><span class="keyword">struct</span> <span class="identifier">push_back_state_and_time</span> 146<span class="special">{</span> 147 <span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span> <span class="identifier">state_type</span> <span class="special">>&</span> <span class="identifier">m_states</span><span class="special">;</span> 148 <span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span> <span class="keyword">double</span> <span class="special">>&</span> <span class="identifier">m_times</span><span class="special">;</span> 149 150 <span class="identifier">push_back_state_and_time</span><span class="special">(</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span> <span class="identifier">state_type</span> <span class="special">></span> <span class="special">&</span><span class="identifier">states</span> <span class="special">,</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">vector</span><span class="special"><</span> <span class="keyword">double</span> <span class="special">></span> <span class="special">&</span><span class="identifier">times</span> <span class="special">)</span> 151 <span class="special">:</span> <span class="identifier">m_states</span><span class="special">(</span> <span class="identifier">states</span> <span class="special">)</span> <span class="special">,</span> <span class="identifier">m_times</span><span class="special">(</span> <span class="identifier">times</span> <span class="special">)</span> <span class="special">{</span> <span class="special">}</span> 152 153 <span class="keyword">void</span> <span class="keyword">operator</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="keyword">double</span> <span class="identifier">t</span> <span class="special">)</span> 154 <span class="special">{</span> 155 <span class="identifier">m_states</span><span class="special">.</span><span class="identifier">push_back</span><span class="special">(</span> <span class="identifier">x</span> <span class="special">);</span> 156 <span class="identifier">m_times</span><span class="special">.</span><span class="identifier">push_back</span><span class="special">(</span> <span class="identifier">t</span> <span class="special">);</span> 157 <span class="special">}</span> 158<span class="special">};</span> 159</pre> 160<p> 161 </p> 162<p> 163 which stores the intermediate steps in a container. Note, the argument structure 164 of the ()-operator: odeint calls the observer exactly in this way, providing 165 the current state and time. Now, you only have to pass this container to 166 the integration function: 167 </p> 168<p> 169</p> 170<pre class="programlisting"><span class="identifier">vector</span><span class="special"><</span><span class="identifier">state_type</span><span class="special">></span> <span class="identifier">x_vec</span><span class="special">;</span> 171<span class="identifier">vector</span><span class="special"><</span><span class="keyword">double</span><span class="special">></span> <span class="identifier">times</span><span class="special">;</span> 172 173<span class="identifier">steps</span> <span class="special">=</span> <span class="identifier">integrate</span><span class="special">(</span> <span class="identifier">harmonic_oscillator</span> <span class="special">,</span> 174 <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="number">0.1</span> <span class="special">,</span> 175 <span class="identifier">push_back_state_and_time</span><span class="special">(</span> <span class="identifier">x_vec</span> <span class="special">,</span> <span class="identifier">times</span> <span class="special">)</span> <span class="special">);</span> 176 177<span class="comment">/* output */</span> 178<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">steps</span><span class="special">;</span> <span class="identifier">i</span><span class="special">++</span> <span class="special">)</span> 179<span class="special">{</span> 180 <span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">times</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special"><<</span> <span class="char">'\t'</span> <span class="special"><<</span> <span class="identifier">x_vec</span><span class="special">[</span><span class="identifier">i</span><span class="special">][</span><span class="number">0</span><span class="special">]</span> <span class="special"><<</span> <span class="char">'\t'</span> <span class="special"><<</span> <span class="identifier">x_vec</span><span class="special">[</span><span class="identifier">i</span><span class="special">][</span><span class="number">1</span><span class="special">]</span> <span class="special"><<</span> <span class="char">'\n'</span><span class="special">;</span> 181<span class="special">}</span> 182</pre> 183<p> 184 </p> 185<p> 186 That is all. You can use functional libraries like <a href="http://www.boost.org/doc/libs/release/libs/lambda/" target="_top">Boost.Lambda</a> 187 or <a href="http://www.boost.org/doc/libs/release/libs/phoenix/" target="_top">Boost.Phoenix</a> 188 to ease the creation of observer functions. 189 </p> 190<p> 191 The full cpp file for this example can be found here: <a href="https://github.com/headmyshoulder/odeint-v2/blob/master/examples/harmonic_oscillator.cpp" target="_top">harmonic_oscillator.cpp</a> 192 </p> 193</div> 194<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr> 195<td align="left"></td> 196<td align="right"><div class="copyright-footer">Copyright © 2009-2015 Karsten Ahnert and Mario Mulansky<p> 197 Distributed under the Boost Software License, Version 1.0. (See accompanying 198 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>) 199 </p> 200</div></td> 201</tr></table> 202<hr> 203<div class="spirit-nav"> 204<a accesskey="p" href="usage__compilation__headers.html"><img src="../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../getting_started.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="../tutorial.html"><img src="../../../../../../../doc/src/images/next.png" alt="Next"></a> 205</div> 206</body> 207</html> 208
