<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
<title>Short Example</title>
<link rel="stylesheet" href="../../boostbook.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.73.2">
<link rel="start" href="../../index.html" title="Chapter 1. boost.sandbox.numeric.odeint">
<link rel="up" href="../getting_started.html" title="Getting started">
<link rel="prev" href="usage__compilation__headers.html" title="Usage, Compilation, Headers">
<link rel="next" href="../tutorial.html" title="Tutorial">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr><td valign="top"></td></tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="usage__compilation__headers.html"><img src="../../images/prev.png" alt="Prev"></a><a accesskey="u" href="../getting_started.html"><img src="../../images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../images/home.png" alt="Home"></a><a accesskey="n" href="../tutorial.html"><img src="../../images/next.png" alt="Next"></a>
</div>
<div class="section" lang="en">
<div class="titlepage"><div><div><h3 class="title">
<a name="boost_sandbox_numeric_odeint.getting_started.short_example"></a><a class="link" href="short_example.html" title="Short Example">Short
Example</a>
</h3></div></div></div>
<p>
Imaging, you want to numerically integrate a harmonic oscillator with friction.
The equations of motion are given by <span class="emphasis"><em>x'' = -x + γ x'</em></span>.
This can be transformed to a system of two coupled first-order differential
equations with new variables <span class="emphasis"><em>x</em></span> and <span class="emphasis"><em>p=x'</em></span>.
To apply numerical integration one first has to design the right hand side
of the equation <span class="emphasis"><em>w' = f(w)</em></span> where in this case <span class="emphasis"><em>w
= (x,p)</em></span>:
</p>
<p>
</p>
<pre class="programlisting"><span class="comment">/* The type of container used to hold the state vector */</span>
<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>
<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>
<span class="comment">/* The rhs of x' = f(x) */</span>
<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="identifier">t</span> <span class="special">)</span>
<span class="special">{</span>
<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>
<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>
<span class="special">}</span>
</pre>
<p>
</p>
<p>
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>
as the state type, but others are also possible, for example <code class="computeroutput"><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="number">2</span><span class="special">></span></code>.
Odeint is designed in such a way that you can easily use your own state types.
Next, we define the ODE which is in this case a simple function. The parameter
signature of this function is crucial: the integration methods 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> <span class="identifier">dxdt</span><span class="special">,</span> <span class="identifier">t</span><span class="special">)</span></code>.
So, even if there is no explicit time dependence, one has to define <code class="computeroutput"><span class="identifier">t</span></code> as a function parameter.
</p>
<p>
Now, we have to define the initial state from which the integration should
start:
</p>
<p>
</p>
<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>
<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><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>
</pre>
<p>
</p>
<p>
For the integration itself we'll use the <code class="computeroutput">integrate</code>
function, which is a convenient way to get quick results. It is based on
the error-controlled <code class="computeroutput">runge_kutta_rk5_ck</code>
stepper (5th order) and uses adaptive stepsize.
</p>
<p>
</p>
<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>
<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>
</pre>
<p>
</p>
<p>
The integrate function expects as parameters the rhs of the ode as defined
above, the initial state <code class="computeroutput"><span class="identifier">x</span></code>,
the start-and end-time of the integration as well as the initial time step.
Note, that <code class="computeroutput">integrate</code>
uses an adaptive stepsize during the integration steps so the time points
will not be equally spaced. The integration returns the number of steps that
were applied.
</p>
<p>
It is, of course, also possible to implement the ode system as a class. The
rhs must then be implemented as a functor having defined the ()-operator:
</p>
<p>
</p>
<pre class="programlisting"><span class="comment">/* The rhs of x' = f(x) defined as a class */</span>
<span class="keyword">class</span> <span class="identifier">harm_osc</span> <span class="special">{</span>
<span class="keyword">double</span> <span class="identifier">m_gam</span><span class="special">;</span>
<span class="keyword">public</span><span class="special">:</span>
<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>
<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="identifier">t</span> <span class="special">)</span>
<span class="special">{</span>
<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>
<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>
<span class="special">}</span>
<span class="special">};</span>
</pre>
<p>
</p>
<p>
which can be used via
</p>
<p>
</p>
<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>
<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>
<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>
</pre>
<p>
</p>
<p>
You surely have already noticed that during the integration a lot of steps
have been done. You might wonder if you could access them do observe the
solution during the iteration. Yes, you can do that. All you have to do is
to provide a reasonable observer. An example is
</p>
<p>
</p>
<pre class="programlisting"><span class="keyword">struct</span> <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="identifier">m_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="identifier">m_times</span><span class="special">;</span>
<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>
<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>
<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>
<span class="special">{</span>
<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>
<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>
<span class="special">}</span>
<span class="special">};</span>
</pre>
<p>
</p>
<p>
which stores the intermediate steps in a container. Now, you only have to
pass this container to the integration function:
</p>
<p>
</p>
<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>
<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>
<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>
<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>
<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>
<span class="comment">/* output */</span>
<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>
<span class="special">{</span>
<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>
<span class="special">}</span>
</pre>
<p>
</p>
<p>
That is all. Of course, you can use functional libraries like <a href="http://www.boost.org/doc/libs/release/doc/html/lambda.html" target="_top">Boost.Lambda</a>
or <a href="http://www.boost.org/doc/libs/1_46_1/libs/spirit/phoenix/doc/html/index.html" target="_top">Boost.Phoenix</a>
to ease the creation of observer functions.
</p>
<p>
The full cpp file for this example can be found here: <a href="../../../../examples/harmonic_oscillator.cpp" target="_top">../../examples/harmonic_oscillator.cpp</a>
</p>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>
<td align="right"><div class="copyright-footer">Copyright © 2009 -2011 Karsten Ahnert and Mario Mulansky<p>
Distributed under the Boost Software License, Version 1.0. (See accompanying
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>)
</p>
</div></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="usage__compilation__headers.html"><img src="../../images/prev.png" alt="Prev"></a><a accesskey="u" href="../getting_started.html"><img src="../../images/up.png" alt="Up"></a><a accesskey="h" href="../../index.html"><img src="../../images/home.png" alt="Home"></a><a accesskey="n" href="../tutorial.html"><img src="../../images/next.png" alt="Next"></a>
</div>
</body>
</html>
Submit an article or tip