Boost-Commit :

Date view	Thread view	Subject view	Author view

Subject: [Boost-commit] svn:boost r59166 - in sandbox/odeint: boost/numeric/odeint libs/numeric/odeint/doc libs/numeric/odeint/doc/html libs/numeric/odeint/doc/html/odeint libs/numeric/odeint/examples
From: mario.mulansky_at_[hidden]
Date: 2010-01-20 15:12:50

Author: mariomulansky
Date: 2010-01-20 15:12:48 EST (Wed, 20 Jan 2010)
New Revision: 59166
URL: http://svn.boost.org/trac/boost/changeset/59166

Log:
started tutorial
some bugfixes according to integrate function
Added:
   sandbox/odeint/libs/numeric/odeint/doc/html/odeint/short_example.html (contents, props changed)
   sandbox/odeint/libs/numeric/odeint/doc/html/odeint/tutorial.html (contents, props changed)
   sandbox/odeint/libs/numeric/odeint/examples/doc_integrate_const.cpp (contents, props changed)
Text files modified:
   sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp | 6 +-
   sandbox/odeint/libs/numeric/odeint/doc/html/index.html | 23 ++++++--
   sandbox/odeint/libs/numeric/odeint/doc/html/odeint/quick_guide.html | 29 +++++++++++
   sandbox/odeint/libs/numeric/odeint/doc/html/standalone_HTML.manifest | 3
   sandbox/odeint/libs/numeric/odeint/doc/odeint.qbk | 97 +++++++++++++++++++++++++++++++++++++++
   sandbox/odeint/libs/numeric/odeint/examples/Jamfile | 3
   sandbox/odeint/libs/numeric/odeint/examples/doc_harm_osc.cpp | 25 ++++++++++
   sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp | 32 +++++++-----
   sandbox/odeint/libs/numeric/odeint/examples/pendulum_vibrating_pivot.cpp | 17 +++---
   9 files changed, 197 insertions(+), 38 deletions(-)

Modified: sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_euler.hpp 2010-01-20 15:12:48 EST (Wed, 20 Jan 2010)
@@ -62,8 +62,8 @@

 template< class DynamicalSystem >
 void do_step( DynamicalSystem &system ,
- container_type &x ,
- const container_type &dxdt ,
+ container_type &x ,
+ const container_type &dxdt ,
 time_type t ,
 time_type dt )
 {
@@ -76,7 +76,7 @@

 template< class DynamicalSystem >
 void do_step( DynamicalSystem &system ,
- container_type &x ,
+ container_type &x ,
 time_type t ,
 time_type dt )
 {

Modified: sandbox/odeint/libs/numeric/odeint/doc/html/index.html
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/doc/html/index.html (original)
+++ sandbox/odeint/libs/numeric/odeint/doc/html/index.html 2010-01-20 15:12:48 EST (Wed, 20 Jan 2010)
@@ -5,7 +5,7 @@
<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="next" href="odeint/quick_guide.html" title="Quick Guide">
+<link rel="next" href="odeint/short_example.html" title="Short Example">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
@@ -17,7 +17,7 @@
<td align="center">More</td>
</tr></table>
<hr>
-<div class="spirit-nav"><a accesskey="n" href="odeint/quick_guide.html"><img src="../../doc/html/images/next.png" alt="Next"></a></div>
+<div class="spirit-nav"><a accesskey="n" href="odeint/short_example.html"><img src="../../doc/html/images/next.png" alt="Next"></a></div>
<div class="chapter" lang="en">
<div class="titlepage"><div>
<div><h2 class="title">
@@ -30,7 +30,7 @@
</h3></div></div>
<div>Copyright © 2009 Karsten Ahnert and Mario Mulansky</div>
<div><div class="legalnotice">
-<a name="id357869"></a>
+<a name="id325363"></a>
 Distributed under the Boost Software License, Version 1.0. (See accompanying
 file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
 
@@ -40,7 +40,16 @@
Table of Contents
<dl>
<dt>Overview</dt>
-<dt>Quick Guide</dt>
+<dt>Short Example</dt>
+<dt>Tutorial</dt>
+<dd><dl>
+<dt>Define the ODE</dt>
+<dt>Stepper Types</dt>
+<dt><a href="odeint/tutorial.html#odeint.tutorial.integration_with_constant_stepsize">Integration
+ with Constant Stepsize</a></dt>
+<dt><a href="odeint/tutorial.html#odeint.tutorial.integration_with_adaptive_stepsize">Integration
+ with Adaptive Stepsize</a></dt>
+</dl></dd>
</dl>
</div>
<div class="section" lang="en">
@@ -61,7 +70,7 @@
 for real life problems and serves just as illustrative example. In odeint, the following algorithms are implemented:
 
<div class="table">
-<a name="id326900"></a>Table 1.1. Stepper Algorithms
+<a name="id326850"></a>Table 1.1. Stepper Algorithms
<div class="table-contents"><table class="table" summary="Stepper Algorithms">
<colgroup>
<col>
@@ -231,10 +240,10 @@
</div>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
-<td align="left">Last revised: January 12, 2010 at 14:09:33 GMT</td>
+<td align="left">Last revised: January 20, 2010 at 19:11:53 GMT</td>
<td align="right"><div class="copyright-footer"></div></td>
</tr></table>
<hr>
-<div class="spirit-nav"><a accesskey="n" href="odeint/quick_guide.html"><img src="../../doc/html/images/next.png" alt="Next"></a></div>
+<div class="spirit-nav"><a accesskey="n" href="odeint/short_example.html"><img src="../../doc/html/images/next.png" alt="Next"></a></div>
</body>
</html>

Modified: sandbox/odeint/libs/numeric/odeint/doc/html/odeint/quick_guide.html
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/doc/html/odeint/quick_guide.html (original)
+++ sandbox/odeint/libs/numeric/odeint/doc/html/odeint/quick_guide.html 2010-01-20 15:12:48 EST (Wed, 20 Jan 2010)
@@ -60,7 +60,7 @@
 state type, but others are also possible, for example tr1/array<double,2>.
 Odeint is designed in such a way that it works with basically any container
 that can be accessed via iterators. The parameter structure of the function
- is crucial: the integration methods will always call them iin the form f(x, dxdt, t). So even if there is no explicit time
+ is crucial: the integration methods will always call them in the form f(x, dxdt, t). So even if there is no explicit time
 dependence, one has to define t as a function
 parameter.
 
@@ -115,6 +115,33 @@
 applied. Note, that in total steps+1 elements will be inserted into times and
 x_t_vec as the initial time and state are inserted before the first step.
 
+
+ It is, of course, also possible to implement the ode system as a class. The
+ rhs must then be defined as the ()-operator:
+ 
+
+ 
+
+
+
+<pre class="programlisting">class harm_osc {
+
+ double m_gam;
+
+public:
+ harm_osc( double gam ) : m_gam(gam) { }
+
+ void operator() (const state_type &x, state_type &dxdt, const double t)
+ {
+ dxdt[0] = x[1];
+ dxdt[1] = -x[0] - m_gam*x[1];
+ }
+};
+</pre>
+
+ 
+
+ 
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>

Added: sandbox/odeint/libs/numeric/odeint/doc/html/odeint/short_example.html
==============================================================================
--- (empty file)
+++ sandbox/odeint/libs/numeric/odeint/doc/html/odeint/short_example.html 2010-01-20 15:12:48 EST (Wed, 20 Jan 2010)
@@ -0,0 +1,160 @@
+<html>
+<head>
+<meta http-equiv="Content-Type" content="text/html; charset=ISO-8859-1">
+<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="../index.html" title="Chapter 1. boost.sandbox.numeric.odeint">
+<link rel="prev" href="../index.html" title="Chapter 1. boost.sandbox.numeric.odeint">
+<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"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../boost.png"></td>
+<td align="center">Home</td>
+<td align="center">Libraries</td>
+<td align="center">People</td>
+<td align="center">FAQ</td>
+<td align="center">More</td>
+</tr></table>
+<hr>
+<div class="spirit-nav">
+<a accesskey="p" href="../index.html"><img src="../../../doc/html/images/prev.png" alt="Prev"></a><a accesskey="u" href="../index.html"><img src="../../../doc/html/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../doc/html/images/home.png" alt="Home"></a><a accesskey="n" href="tutorial.html"><img src="../../../doc/html/images/next.png" alt="Next"></a>
+</div>
+<div class="section" lang="en">
+<div class="titlepage"><div><div><h2 class="title" style="clear: both">
+<a name="odeint.short_example"></a><a class="link" href="short_example.html" title="Short Example">Short Example</a>
+</h2></div></div></div>
+
+ This section gives a quick introduction to the most important features of the
+ library using a number of instructive examples. Image, for example, you want
+ to numerically integrate a harmonic oscillator with friction. The equations
+ of motion are given by x'' = -x + gamma x'. This can be transformed to a system
+ of two first-order differential equations with new variables x and p=x'. To
+ apply numerical integration one first has to design the right hand side of
+ the equation w' = f(w) where in this case w = (x,p):
+ 
+
+ 
+
+
+
+<pre class="programlisting">/* The type of container used to hold the state vector */
+typedef std::vector<double> state_type;
+
+const double gam = 0.15;
+
+/* The rhs of x' = f(x) */
+void harmonic_oscillator(const state_type &x, state_type &dxdt, const double t)
+{
+ dxdt[0] = x[1];
+ dxdt[1] = -x[0] - gam*x[1];
+}
+</pre>
+
+ 
+
+ 
+
+ Here we chose vector<double> as the
+ state type, but others are also possible, for example tr1/array<double,2>.
+ Odeint is designed in such a way that it works with basically any container
+ that can be accessed via iterators. The parameter structure of the function
+ is crucial: the integration methods will always call them in the form f(x, dxdt, t). So even if there is no explicit time
+ dependence, one has to define t as a function
+ parameter.
+ 
+
+ Now, we have to define the initial state from which the integration should
+ start:
+ 
+
+ 
+
+
+
+<pre class="programlisting">state_type x(2);
+x[0] = 1.0; // start at x=1.0, p=0.0
+x[1] = 0.0;
+</pre>
+
+ 
+
+ 
+
+ 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. The results are stored into
+ two InsertIterators (time and state) that must be provided.
+ 
+
+ 
+
+
+
+<pre class="programlisting">vector<double> times;
+vector<state_type> x_t_vec;
+
+size_t steps = integrate( harmonic_oscillator ,
+ x , 0.0 , 10.0 ,
+ back_inserter( times ) ,
+ back_inserter( x_t_vec ) );
+</pre>
+
+ 
+
+ 
+
+ The integrate function expects as parameters the rhs of the ode as defined
+ above, the initial state x, the start- and end-time of the integration and
+ two InsertIterators for the times and states where the current time and the
+ current state after each timestep is stored. 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. Note, that in total steps+1 elements will be inserted into times and
+ x_t_vec as the initial time and state are inserted before the first step.
+ 
+
+ It is, of course, also possible to implement the ode system as a class. The
+ rhs must then be defined as the ()-operator:
+ 
+
+ 
+
+
+
+<pre class="programlisting">class harm_osc {
+
+ double m_gam;
+
+public:
+ harm_osc( double gam ) : m_gam(gam) { }
+
+ void operator() (const state_type &x, state_type &dxdt, const double t)
+ {
+ dxdt[0] = x[1];
+ dxdt[1] = -x[0] - m_gam*x[1];
+ }
+};
+</pre>
+
+ 
+
+ 
+</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 Karsten Ahnert and Mario Mulansky
+ Distributed under the Boost Software License, Version 1.0. (See accompanying
+ file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ 
+</div></td>
+</tr></table>
+<hr>
+<div class="spirit-nav">
+<a accesskey="p" href="../index.html"><img src="../../../doc/html/images/prev.png" alt="Prev"></a><a accesskey="u" href="../index.html"><img src="../../../doc/html/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../doc/html/images/home.png" alt="Home"></a><a accesskey="n" href="tutorial.html"><img src="../../../doc/html/images/next.png" alt="Next"></a>
+</div>
+</body>
+</html>

Added: sandbox/odeint/libs/numeric/odeint/doc/html/odeint/tutorial.html
==============================================================================
--- (empty file)
+++ sandbox/odeint/libs/numeric/odeint/doc/html/odeint/tutorial.html 2010-01-20 15:12:48 EST (Wed, 20 Jan 2010)
@@ -0,0 +1,367 @@
+<html>
+<head>
+<meta http-equiv="Content-Type" content="text/html; charset=ISO-8859-1">
+<title>Tutorial</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="../index.html" title="Chapter 1. boost.sandbox.numeric.odeint">
+<link rel="prev" href="short_example.html" title="Short Example">
+</head>
+<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
+<table cellpadding="2" width="100%"><tr>
+<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../boost.png"></td>
+<td align="center">Home</td>
+<td align="center">Libraries</td>
+<td align="center">People</td>
+<td align="center">FAQ</td>
+<td align="center">More</td>
+</tr></table>
+<hr>
+<div class="spirit-nav">
+<a accesskey="p" href="short_example.html"><img src="../../../doc/html/images/prev.png" alt="Prev"></a><a accesskey="u" href="../index.html"><img src="../../../doc/html/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../doc/html/images/home.png" alt="Home"></a>
+</div>
+<div class="section" lang="en">
+<div class="titlepage"><div><div><h2 class="title" style="clear: both">
+<a name="odeint.tutorial"></a><a class="link" href="tutorial.html" title="Tutorial">Tutorial</a>
+</h2></div></div></div>
+<div class="toc"><dl>
+<dt>Define the ODE</dt>
+<dt>Stepper Types</dt>
+<dt><a href="tutorial.html#odeint.tutorial.integration_with_constant_stepsize">Integration
+ with Constant Stepsize</a></dt>
+<dt><a href="tutorial.html#odeint.tutorial.integration_with_adaptive_stepsize">Integration
+ with Adaptive Stepsize</a></dt>
+</dl></div>
+<div class="section" lang="en">
+<div class="titlepage"><div><div><h3 class="title">
+<a name="odeint.tutorial.define_the_ode"></a><a class="link" href="tutorial.html#odeint.tutorial.define_the_ode" title="Define the ODE">Define the ODE</a>
+</h3></div></div></div>
+
+ First of all, you have to specify the datatype that represents a state of
+ your system x. Mathematically, this usually
+ is an n-dimensional vector with real numbers or complex numbers as scalar
+ objects. For odeint the most natural way is to use <code class="computeroutput">vector< double ></code> or <code class="computeroutput">vector< complex< double > ></code>
+ to represent the system state. However, odeint can deal with other container
+ types than vector as well, e.g. <code class="computeroutput">tr1/array< double , N ></code>
+ as long as it is able to obtain a ForwardIterator going through all of the
+ container's elements. The scalar type must have several operators ( +, -,
+ +<code class="literal">, -</code> ) and the <code class="computeroutput">abs()</code>-fcuntion defined . Furthermore, one can
+ choose the datatype of the time (that is, the parameter to which respect
+ the differentiation is done). The standard type for this is <code class="computeroutput">double</code>, but one might be able to choose, for
+ example, <code class="computeroutput">complex<
+ double ></code>
+ as well (untested). It must be possible to multiply the time type and the
+ scalar type of the vector. For the most cases <code class="computeroutput">vector< double ></code> as state type and the standard <code class="computeroutput">double</code> for the time type should be sufficient.
+ 
+
+ To integrate a differential equation numerically, one has to define the rhs
+ of the equation x' = f(x). In odeint you supply this
+ function in terms of an object that implements the ()-operator with a certain
+ parameter structure. Hence, the straight forward way would be to just define
+ a function, e.g:
+ 
+
+
+
+<pre class="programlisting">/* The type of container used to hold the state vector */
+typedef std::vector<double> state_type;
+
+const double gam = 0.15;
+
+/* The rhs of x' = f(x) */
+void harmonic_oscillator(const state_type &x, state_type &dxdt, const double t)
+{
+ dxdt[0] = x[1];
+ dxdt[1] = -x[0] - gam*x[1];
+}
+</pre>
+
+ 
+
+ The parameters of the function must follow the example above where x is the current state, dxdt
+ is the derivative x' and should be filled by the function
+ with f(x) and t is
+ the current time.
+ 
+
+ A more sophisticated approach is to implement the system as a class where
+ the rhs function is defined as the ()-operator of the class with the same
+ parameter structure as above:
+ 
+
+
+
+<pre class="programlisting">class harm_osc {
+
+ double m_gam;
+
+public:
+ harm_osc( double gam ) : m_gam(gam) { }
+
+ void operator() (const state_type &x, state_type &dxdt, const double t)
+ {
+ dxdt[0] = x[1];
+ dxdt[1] = -x[0] - m_gam*x[1];
+ }
+};
+</pre>
+
+ 
+
+ odeint can deal with instances of such classes instead of pure functions
+ which allows for cleaner code.
+ 
+</div>
+<div class="section" lang="en">
+<div class="titlepage"><div><div><h3 class="title">
+<a name="odeint.tutorial.stepper_types"></a><a class="link" href="tutorial.html#odeint.tutorial.stepper_types" title="Stepper Types">Stepper Types</a>
+</h3></div></div></div>
+
+ Numerical integration works iteratevly, that means you start at a state
+ x(t) and performs a timestep of length dt
+ to obtain the approximate state x(t+dt). There exist
+ many different methods to perform such a timestep each of which has a certain
+ order q. If the order of a method is q
+ than it is accurate up to term ~dt^q that means the
+ error in x made by such a step is ~dt^(q+1).
+ odeint provides several steppers of different orders from which you can choose.
+ There are three types of steppers: Stepper,
+ ErrorStepper and ControlledStepper.
+ 
+<div class="table">
+<a name="id368578"></a>Table 1.2. Stepper Algorithms
+<div class="table-contents"><table class="table" summary="Stepper Algorithms">
+<colgroup>
+<col>
+<col>
+<col>
+<col>
+</colgroup>
+<thead><tr>
+<th>
+ 
+ Method
+ 
+ </th>
+<th>
+ 
+ Class
+ 
+ </th>
+<th>
+ 
+ Order
+ 
+ </th>
+<th>
+ 
+ Error Estimation
+ 
+ </th>
+</tr></thead>
+<tbody>
+<tr>
+<td>
+ 
+ Euler
+ 
+ </td>
+<td>
+ 
+ stepper_euler
+ 
+ </td>
+<td>
+ 
+ 1
+ 
+ </td>
+<td>
+ 
+ No
+ 
+ </td>
+</tr>
+<tr>
+<td>
+ 
+ Runge-Kutta 4
+ 
+ </td>
+<td>
+ 
+ stepper_rk4
+ 
+ </td>
+<td>
+ 
+ 4
+ 
+ </td>
+<td>
+ 
+ No
+ 
+ </td>
+</tr>
+<tr>
+<td>
+ 
+ Runge-Kutta Cash-Karp
+ 
+ </td>
+<td>
+ 
+ stepper_rk5_ck
+ 
+ </td>
+<td>
+ 
+ 5
+ 
+ </td>
+<td>
+ 
+ Yes (Order 4)
+ 
+ </td>
+</tr>
+<tr>
+<td>
+ 
+ Runge-Kutta Fehlberg
+ 
+ </td>
+<td>
+ 
+ stepper_rk78_fehlberg
+ 
+ </td>
+<td>
+ 
+ 7
+ 
+ </td>
+<td>
+ 
+ Yes (Order 8)
+ 
+ </td>
+</tr>
+<tr>
+<td>
+ 
+ Midpoint
+ 
+ </td>
+<td>
+ 
+ stepper_midpoint
+ 
+ </td>
+<td>
+ 
+ variable
+ 
+ </td>
+<td>
+ 
+ No
+ 
+ </td>
+</tr>
+<tr>
+<td>
+ 
+ Bulirsch-Stoer
+ 
+ </td>
+<td>
+ 
+ controlled_stepper_bs
+ 
+ </td>
+<td>
+ 
+ variable
+ 
+ </td>
+<td>
+ 
+ Controlled
+ 
+ </td>
+</tr>
+</tbody>
+</table></div>
+</div>
+ 
+</div>
+<div class="section" lang="en">
+<div class="titlepage"><div><div><h3 class="title">
+<a name="odeint.tutorial.integration_with_constant_stepsize"></a><a class="link" href="tutorial.html#odeint.tutorial.integration_with_constant_stepsize" title="Integration with Constant Stepsize">Integration
+ with Constant Stepsize</a>
+</h3></div></div></div>
+
+ The basic stepper just performs one timestep and doesn't give you any information
+ about the error that was made (except that you know it is of order q). Such
+ steppers are used with constant stepsize that should be chosen small enough
+ to have reasonable small errors. However, you should apply some sort of validity
+ check of your results (such as observing conserved quantities) becasue you
+ have no other control of the error. The following example defines basic stepper
+ based on the classical Runge-Kutta scheme of 4th order.
+ 
+
+ 
+
+
+
+<pre class="programlisting">stepper_rk4< state_type > rk4;
+</pre>
+
+ 
+
+ 
+
+ The declaration of the stepper requires the state type as template parameter.
+ The integration can now be done by using the <code class="computeroutput">integrate_const( Stepper, System, start_time, step_size, state, end_time)</code> function from odeint:
+ 
+
+ 
+
+
+
+<pre class="programlisting">integrate_const( rk4, harmonic_oscillator, 0.0, 0.01, x, 10.0 );
+</pre>
+
+ 
+
+ 
+
+ This call integrates the system defined by <code class="computeroutput">harmonic_oscillator</code>
+ using the rk4 method from t=0 to 10 with a stepsize dt=0.01 with the initial
+ condition given in <code class="computeroutput">x</code>. The
+ result, x(t=10) is stored in <code class="computeroutput">x</code>
+ (in-place).
+ 
+</div>
+<div class="section" lang="en"><div class="titlepage"><div><div><h3 class="title">
+<a name="odeint.tutorial.integration_with_adaptive_stepsize"></a><a class="link" href="tutorial.html#odeint.tutorial.integration_with_adaptive_stepsize" title="Integration with Adaptive Stepsize">Integration
+ with Adaptive Stepsize</a>
+</h3></div></div></div></div>
+</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 Karsten Ahnert and Mario Mulansky
+ Distributed under the Boost Software License, Version 1.0. (See accompanying
+ file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ 
+</div></td>
+</tr></table>
+<hr>
+<div class="spirit-nav">
+<a accesskey="p" href="short_example.html"><img src="../../../doc/html/images/prev.png" alt="Prev"></a><a accesskey="u" href="../index.html"><img src="../../../doc/html/images/up.png" alt="Up"></a><a accesskey="h" href="../index.html"><img src="../../../doc/html/images/home.png" alt="Home"></a>
+</div>
+</body>
+</html>

Modified: sandbox/odeint/libs/numeric/odeint/doc/html/standalone_HTML.manifest
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/doc/html/standalone_HTML.manifest (original)
+++ sandbox/odeint/libs/numeric/odeint/doc/html/standalone_HTML.manifest 2010-01-20 15:12:48 EST (Wed, 20 Jan 2010)
@@ -1,2 +1,3 @@
index.html
-odeint/quick_guide.html
+odeint/short_example.html
+odeint/tutorial.html

Modified: sandbox/odeint/libs/numeric/odeint/doc/odeint.qbk
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/doc/odeint.qbk (original)
+++ sandbox/odeint/libs/numeric/odeint/doc/odeint.qbk 2010-01-20 15:12:48 EST (Wed, 20 Jan 2010)
@@ -43,7 +43,7 @@

[endsect]

-[section Quick Guide]
+[section Short Example]
This section gives a quick introduction to the most important features of the
library using a number of instructive examples.
Image, for example, you want to numerically integrate a harmonic oscillator
@@ -60,7 +60,7 @@
for example [*tr1/array<double,2>]. Odeint is designed in such a way that it
works with basically any container that can be accessed via iterators.
The parameter structure of the function is crucial: the integration methods will
-always call them iin the form [*f(x, dxdt, t)]. So even if there is no explicit
+always call them in the form [*f(x, dxdt, t)]. So even if there is no explicit
time dependence, one has to define [*t] as a function parameter.

Now, we have to define the initial state from which the integration should start:
@@ -88,4 +88,97 @@
Note, that in total steps+1 elements will be inserted into times and x_t_vec as
the initial time and state are inserted before the first step.

+It is, of course, also possible to implement the ode system as a class. The rhs
+must then be defined as the ()-operator:
+
+[rhs_class]
+
+[endsect]
+
+[section Tutorial]
+
+[section Define the ODE]
+First of all, you have to specify the datatype that represents a state of your
+system [*x]. Mathematically, this usually is an n-dimensional vector with
+real numbers or complex numbers as scalar objects. For odeint the most natural
+way is to use `vector< double >` or `vector< complex< double > >` to represent
+the system state. However, odeint can deal with other container types than
+vector as well, e.g. `tr1/array< double , N >` as long as it is able to obtain
+a ForwardIterator going through all of the container's elements. The scalar type
+must have several operators ( +, -, +=, -= ) and the `abs()`-fcuntion defined .
+Furthermore, one can choose the datatype of the time (that is, the parameter to
+which respect the differentiation is done). The standard type for this is
+`double`, but one might be able to choose, for example, `complex< double >` as
+well (untested). It must be possible to multiply the time type and the scalar
+type of the vector. For the most cases `vector< double >` as state type and the
+standard `double` for the time type should be sufficient.
+
+To integrate a differential equation numerically, one has to define the rhs of
+the equation ['x' = f(x)]. In odeint you supply this function in terms of
+an object that implements the ()-operator with a certain parameter structure.
+Hence, the straight forward way would be to just define a function, e.g:
+[rhs_function]
+The parameters of the function must follow the example above where [*x] is the
+current state, [*dxdt] is the derivative ['x'] and should be filled by the
+function with ['f(x)] and [*t] is the current time.
+
+A more sophisticated approach is to implement the system as a class where the
+rhs function is defined as the ()-operator of the class with the same parameter
+structure as above:
+[rhs_class]
+odeint can deal with instances of such classes instead of pure functions which
+allows for cleaner code.
+
+[endsect]
+
+[section Stepper Types]
+Numerical integration works iteratevly, that means you start at a state ['x(t)]
+and performs a timestep of length ['dt] to obtain the approximate state
+['x(t+dt)]. There exist many different methods to perform such a timestep each
+of which has a certain order ['q]. If the order of a method is ['q] than it is
+accurate up to term ['~dt^q] that means the error in ['x] made by such a step
+is ['~dt^(q+1)]. odeint provides several steppers of different orders from which
+you can choose. There are three types of steppers: [*Stepper], [*ErrorStepper]
+and [*ControlledStepper].
+
+[table Stepper Algorithms
+ [[Method] [Class] [Order] [Error Estimation]]
+ [[Euler] [stepper_euler] [1] [No]]
+ [[Runge-Kutta 4] [stepper_rk4] [4] [No]]
+ [[Runge-Kutta Cash-Karp] [stepper_rk5_ck] [5] [Yes (Order 4)]]
+ [[Runge-Kutta Fehlberg] [stepper_rk78_fehlberg] [7] [Yes (Order 8)]]
+ [[Midpoint] [stepper_midpoint] [variable] [No]]
+ [[Bulirsch-Stoer] [controlled_stepper_bs] [variable] [Controlled]]
+]
+
+[endsect]
+
+[section Integration with Constant Stepsize]
+[import ../examples/doc_integrate_const.cpp]
+The basic stepper just performs one timestep and doesn't give you any
+information about the error that was made (except that you know it is of order
+q). Such steppers are used with constant stepsize that should be chosen small
+enough to have reasonable small errors. However, you should apply some sort of
+validity check of your results (such as observing conserved quantities) becasue
+you have no other control of the error. The following example defines basic
+stepper based on the classical Runge-Kutta scheme of 4th order.
+
+[define_stepper]
+
+The declaration of the stepper requires the state type as template parameter.
+The integration can now be done by using the `integrate_const( Stepper, System,
+ start_time, step_size, state, end_time)` function from odeint:
+
+[integrate]
+
+This call integrates the system defined by `harmonic_oscillator` using the rk4
+method from t=0 to 10 with a stepsize dt=0.01 with the initial condition given
+in `x`. The result, ['x(t=10)] is stored in `x` (in-place).
+
+[endsect]
+
+[section Integration with Adaptive Stepsize]
+
+[endsect]
+
[endsect]
\ No newline at end of file

Modified: sandbox/odeint/libs/numeric/odeint/examples/Jamfile
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/Jamfile (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/Jamfile 2010-01-20 15:12:48 EST (Wed, 20 Jan 2010)
@@ -26,4 +26,5 @@
exe pendulum_vibrating_pivot : pendulum_vibrating_pivot.cpp ;
exe dnls_stepper_compare : dnls_stepper_compare.cpp ;

-exe doc_harm_osc : doc_harm_osc.cpp ;
\ No newline at end of file
+exe doc_harm_osc : doc_harm_osc.cpp ;
+exe doc_integrate_const : doc_integrate_const.cpp ;
\ No newline at end of file

Modified: sandbox/odeint/libs/numeric/odeint/examples/doc_harm_osc.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/doc_harm_osc.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/doc_harm_osc.cpp 2010-01-20 15:12:48 EST (Wed, 20 Jan 2010)
@@ -16,6 +16,22 @@
}
//]

+//[ rhs_class
+class harm_osc {
+
+ double m_gam;
+
+public:
+ harm_osc( double gam ) : m_gam(gam) { }
+
+ void operator() (const state_type &x, state_type &dxdt, const double t)
+ {
+ dxdt[0] = x[1];
+ dxdt[1] = -x[0] - m_gam*x[1];
+ }
+};
+//]
+
int main(int argc, char **argv)
{
 using namespace std;
@@ -37,6 +53,15 @@
 back_inserter( x_t_vec ) );
 //]

+ /* the same as above using the class */
+ /*
+ harm_osc ho(0.15);
+ steps = integrate( ho ,
+ x , 0.0 , 10.0 ,
+ back_inserter( times ) ,
+ back_inserter( x_t_vec ) );
+ */
+
 //[ output
 for( size_t i=0; i<=steps; i++ ) { //initial state is 0th entry
 cout << times[i] << '\t' << x_t_vec[i][0] << '\t' << x_t_vec[i][1] << '\n';

Added: sandbox/odeint/libs/numeric/odeint/examples/doc_integrate_const.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/libs/numeric/odeint/examples/doc_integrate_const.cpp 2010-01-20 15:12:48 EST (Wed, 20 Jan 2010)
@@ -0,0 +1,45 @@
+#include <iostream>
+
+#include <boost/numeric/odeint.hpp>
+
+using namespace std;
+
+/* The type of container used to hold the state vector */
+typedef std::vector<double> state_type;
+
+class harm_osc {
+
+ double m_gam;
+
+public:
+ harm_osc( double gam ) : m_gam(gam) { }
+
+ void operator() (const state_type &x, state_type &dxdt, const double t)
+ {
+ dxdt[0] = x[1];
+ dxdt[1] = -x[0] - m_gam*x[1];
+ }
+};
+
+int main(int argc, char **argv)
+{
+ using namespace std;
+ using namespace boost::numeric::odeint;
+
+ state_type x(2);
+ x[0] = 1.0; // start at x=1.0, p=0.0
+ x[1] = 0.0;
+
+ harm_osc harmonic_oscillator(0.15);
+
+ //[ define_stepper
+ stepper_rk4< state_type > rk4;
+ //]
+
+ //[ integrate
+ integrate_const( rk4, harmonic_oscillator, 0.0, 0.01, x, 10.0 );
+ //]
+
+ cout << x[0] << '\t' << x[1] << endl;
+
+}

Modified: sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp 2010-01-20 15:12:48 EST (Wed, 20 Jan 2010)
@@ -60,25 +60,28 @@

 vector<state_type> x1_t_vec;
+ vector<double> t1_vec;
 vector<state_type> x2_t_vec;
+ vector<double> t2_vec;
 vector<state_type> x3_t_vec;
- vector<double> times(time_points);
- for( size_t i=0; i<time_points; i++ ) {
- times[i] = 0.1*i;
- }
+ vector<double> t3_vec;

 stepper_half_step< stepper_euler< state_type > > euler;
 controlled_stepper_standard< stepper_half_step< stepper_euler< state_type > > >
 euler_controlled( euler , eps_abs, eps_rel, 1.0, 1.0);
 size_t steps = integrate( euler_controlled, lorenz, x1,
- times, 1E-4, back_inserter(x1_t_vec));
+ 0.0, 10.0, 1E-4,
+ back_inserter(t1_vec),
+ back_inserter(x1_t_vec));

 clog << "Euler Half Step: #steps " << steps << endl;

 stepper_half_step< stepper_rk4< state_type > > rk4;
 controlled_stepper_standard< stepper_half_step< stepper_rk4< state_type > > >
 rk4_controlled( rk4 , eps_abs, eps_rel, 1.0, 1.0);
- steps = integrate( rk4_controlled, lorenz, x2, times, 1E-4, back_inserter(x2_t_vec));
+ steps = integrate( rk4_controlled, lorenz, x2, 0.0, 10.0, 1E-4,
+ back_inserter(t2_vec),
+ back_inserter(x2_t_vec));

 clog << "RK4 Half Step: #steps " << steps << endl;

@@ -86,21 +89,22 @@
 stepper_rk5_ck< state_type > rk5;
 controlled_stepper_standard< stepper_rk5_ck< state_type > >
 rk5_controlled( rk5 , eps_abs, eps_rel, 1.0, 1.0);
- steps = integrate( rk5_controlled, lorenz, x3, times, 1E-4, back_inserter(x3_t_vec));
+ steps = integrate( rk5_controlled, lorenz, x3, 0.0, 10.0, 1E-4,
+ back_inserter(t3_vec),
+ back_inserter(x3_t_vec));

 clog << "RK5 Cash-Karp: #steps: " << steps << endl;

 cout.precision(16);
 cout.setf(ios::fixed,ios::floatfield);
+
+ cout << t1_vec.size() << tab << t2_vec.size() << tab << t3_vec.size() << endl;


- for( size_t i=0; i<time_points; i++ ) {
- //cout << "current state: " ;
- cout << times[i] << tab;
- cout << x1_t_vec[i][0] << tab << x1_t_vec[i][1] << tab << x1_t_vec[i][2] << tab;
- cout << x2_t_vec[i][0] << tab << x2_t_vec[i][1] << tab << x2_t_vec[i][2] << tab;
- cout << x3_t_vec[i][0] << tab << x3_t_vec[i][1] << tab << x3_t_vec[i][2] << endl;
- }
+ //cout << "current state: " ;
+ cout << (x1_t_vec.back())[0] << tab << (x1_t_vec.back())[1] << tab << (x1_t_vec.back())[2] << tab;
+ cout << x2_t_vec.back()[0] << tab << x2_t_vec.back()[1] << tab << x2_t_vec.back()[2] << tab;
+ cout << x3_t_vec.back()[0] << tab << x3_t_vec.back()[1] << tab << x3_t_vec.back()[2] << endl;

 return 0;
}

Modified: sandbox/odeint/libs/numeric/odeint/examples/pendulum_vibrating_pivot.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/pendulum_vibrating_pivot.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/pendulum_vibrating_pivot.cpp 2010-01-20 15:12:48 EST (Wed, 20 Jan 2010)
@@ -33,9 +33,8 @@
typedef std::tr1::array< double , 2 > state_type;

const double alpha = 0.1;
-const double omega = 10;
+const double omega = 20;
const double a = 0.1;
-const size_t time_points = 1000;

/*
 Defines the right hand side f(x,t) of the dynamical equations dx/dt = f(x,t)
@@ -51,10 +50,7 @@
{
 state_type x = {{ 1.0, 0.0 }};

- vector<double> times(time_points);
- for( size_t i=0; i<time_points; i++ ) {
- times[i] = 0.1*i;
- }
+ vector<double> times;
 vector<state_type> x_t_vec;

 stepper_half_step< stepper_rk4< state_type > > stepper;
@@ -62,7 +58,10 @@
 controlled_stepper_standard< stepper_half_step< stepper_rk4< state_type > > >
 controlled_stepper( stepper, 1E-6 , 1E-7 , 1.0 , 1.0 );

- size_t steps = integrate( controlled_stepper, my_system, x, times, 1E-4, back_inserter(x_t_vec)); \
+ size_t steps = integrate( controlled_stepper, my_system, x,
+ 0.0, 100.0,1E-4,
+ back_inserter(times),
+ back_inserter(x_t_vec));

 clog << "Steps: " << steps << endl;

@@ -70,10 +69,10 @@
 cout.setf(ios::fixed,ios::floatfield);


- for( size_t i=0; i<time_points; i++ ) {
+ for( size_t i=0; i<times.size(); i++ ) {
 //cout << "current state: " ;
 cout << times[i] << tab;
- cout << x_t_vec[i][0] << tab << x_t_vec[i][1] << tab << x_t_vec[i][2] << endl;
+ cout << x_t_vec[i][0] << tab << x_t_vec[i][1] << endl;
 }

 return 0;

Date view	Thread view	Subject view	Author view

Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk