Boost-Commit :

Date view	Thread view	Subject view	Author view

Subject: [Boost-commit] svn:boost r63652 - in sandbox/odeint/libs/numeric/odeint: doc doc/html doc/html/boost_numeric_odeint examples
From: karsten.ahnert_at_[hidden]
Date: 2010-07-05 10:50:12

Next message: karsten.ahnert_at_[hidden]: "[Boost-commit] svn:boost r63653 - in sandbox/odeint/libs/numeric/odeint/doc: . html html/boost_numeric_odeint"
Previous message: anthony_at_[hidden]: "[Boost-commit] svn:boost r63651 - trunk/libs/thread/test"

Author: karsten
Date: 2010-07-05 10:50:10 EDT (Mon, 05 Jul 2010)
New Revision: 63652
URL: http://svn.boost.org/trac/boost/changeset/63652

Log:
changes in the solar system example tutorial
Text files modified:
   sandbox/odeint/libs/numeric/odeint/doc/html/boost_numeric_odeint/concepts.html | 2
   sandbox/odeint/libs/numeric/odeint/doc/html/boost_numeric_odeint/short_example.html | 20 ++--
   sandbox/odeint/libs/numeric/odeint/doc/html/boost_numeric_odeint/tutorial.html | 157 +++++++++++++++++++++++++++++++++------
   sandbox/odeint/libs/numeric/odeint/doc/html/index.html | 9 +-
   sandbox/odeint/libs/numeric/odeint/doc/tutorial.qbk | 41 +++++----
   sandbox/odeint/libs/numeric/odeint/examples/point_type.hpp | 10 ++
   sandbox/odeint/libs/numeric/odeint/examples/solar_system.cpp | 37 +++++---
   7 files changed, 203 insertions(+), 73 deletions(-)

Modified: sandbox/odeint/libs/numeric/odeint/doc/html/boost_numeric_odeint/concepts.html
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/doc/html/boost_numeric_odeint/concepts.html (original)
+++ sandbox/odeint/libs/numeric/odeint/doc/html/boost_numeric_odeint/concepts.html 2010-07-05 10:50:10 EDT (Mon, 05 Jul 2010)
@@ -304,7 +304,7 @@

 Error steppers execute one timestep of a specific order with a given stepsize.
 Additionally, an error estimation of the obtained result is computed that
- can be used to control the error introduced by the time discretization. As
+ can be used to control the error introduced by the time discretization. Like
 the basic steppers, error steppers usually allocate internal memory to store
 intermediate function call results. If state types with variable size are
 used (e.g. <code class="computeroutput">vector</code>), it has

Modified: sandbox/odeint/libs/numeric/odeint/doc/html/boost_numeric_odeint/short_example.html
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/doc/html/boost_numeric_odeint/short_example.html (original)
+++ sandbox/odeint/libs/numeric/odeint/doc/html/boost_numeric_odeint/short_example.html 2010-07-05 10:50:10 EDT (Mon, 05 Jul 2010)
@@ -21,11 +21,11 @@
</h2></div></div></div>

 This section gives a quick introduction to the most important features of the
- library. 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 design
- the right hand side of the equation w' = f(w) where in this case with new variables
- x and p=x'. To apply numerical integration one first has to w = (x,p):
+ library. Imagine, 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
+ by introducing the new variable p = x'. If we now use w = (x,p) the equation
+ of motion above can be written as w' = f(w):
 

 
@@ -52,10 +52,12 @@
 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.
+ that can be accessed via iterators. Moreover, any structure can be used as
+ state type, provided vector-space operations are provided for them. 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

Modified: sandbox/odeint/libs/numeric/odeint/doc/html/boost_numeric_odeint/tutorial.html
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/doc/html/boost_numeric_odeint/tutorial.html (original)
+++ sandbox/odeint/libs/numeric/odeint/doc/html/boost_numeric_odeint/tutorial.html 2010-07-05 10:50:10 EDT (Mon, 05 Jul 2010)
@@ -55,13 +55,12 @@
 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>
+ However, odeint can deal with other container types 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
+ the container's elements. The scalar type must have several operations
+ ( + , - , += , -= ) and the <code class="computeroutput">abs()</code>-function 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,
+ the differentiation is done). The standard type for this is <code class="computeroutput">double</code>, but it should be possible to use,
 for example, <code class="computeroutput">complex<
 double ></code>
 as well (untested). It must be possible to multiply the time type and the
@@ -140,7 +139,7 @@
 Types</a>
</h4></div></div></div>

- Numerical integration works iteratevly, that means you start at a state
+ Numerical integration works iteratively, 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
@@ -380,7 +379,7 @@
 , a_dxdt
 )</code> function that takes an error stepper
 as parameter and four values defining the maximal absolute and relative
- error allowed for on integration step. The standard controlled stepper
+ error allowed for one integration step. The standard controlled stepper
 created by this method ensures that the error err
 of the solution fulfills err < eps_abs + eps_rel * ( a_x *
 |x| + a_dxdt * dt * |dxdt| ) by decreasesing the step size.
@@ -435,7 +434,7 @@
 derived from the Hamiltonian
 

- H = sum over i pi^2 / 2 + sum j V( qi , qj )
+ H = sum over i pi^2 / 2 m_i + sum j V( qi , qj )
 

 via dqi = dH / dpi, dpi = - dH / dq1. V(qi,qj) is the interaction potential.
@@ -443,8 +442,9 @@

 In time independent Hamiltonian system the energy is conserved and special
 integration methods have to be applied in order to ensure energy conservation.
- The odeint library provides two stepper classes for Hamiltonian systems,
- which are separable and can be written in the form H = sum pi^2/2 + Hq.
+ The odeint library provides classes for Hamiltonian systems, which are
+ separable and can be written in the form H = sum pi^2/2 mi + Hq(q), where
+ Hq(q) only depends on the coordinates.
 

 hamiltonian_stepper_euler hamiltonian_stepper_rk
@@ -452,7 +452,7 @@

 Alltough this functional form might look a bit arbitrary it covers nearly
 all classical mechanical systems with inertia and without dissipation,
- or where the equations of motion can be written in the form dqi=pi dpi=f(qi).
+ or where the equations of motion can be written in the form dqi=mi pi dpi=f(qi).
 
</div>
<div class="section" lang="en">
@@ -465,29 +465,109 @@
 space as well as the velocity. Therefore, we use the operators from <boost/operator.hpp>:
 

- show the code
+ 
+
+
+
+<pre class="programlisting">template< class T , size_t Dim >
+class point :
+ boost::additive1< point< T , Dim > ,
+ boost::additive2< point< T , Dim > , T ,
+ boost::multiplicative2< point< T , Dim > , T
+ > > >
+{
+ // ...
+};
+
+// define the scalar_prod(), norm() and abs()
+</pre>
+
+ 
+
 

- The next step is to define the state type and the system (derivative) function.
- As state type we use std::tr1::array and a state type represents all space
- coordinates q or all momenta coordinates p. As system function we have
- to provide f(q)
+ The next step is to define a container type storing the values of q and
+ p and to define systems (derivative) functions. As container type we use
+ std::tr1::array and the state type is then simply
+ 
+
+ 
+
+
+
+<pre class="programlisting">typedef point< double , 3 > point_type;
+typedef std::tr1::array< point_type , n > container_type;
+typedef std::tr1::array< double , n > mass_type;
+typedef pair< container_type , container_type > state_type;
+</pre>
+
+ 
+
 

- show the code
+ and represents all space coordinates q or all momenta coordinates p. As
+ system function we have to provide f(p) = dq and f(q) = -dp, which acts
+ only on p or q:
 

- Note, that we have allready define the masses of all planets in the solar
- system.
+ 
+
+
+
+<pre class="programlisting">struct solar_system_momentum
+{
+ const mass_type &m_masses;
+
+ solar_system_momentum( const mass_type &masses ) : m_masses( masses ) { }
+
+ void operator()( const container_type &p , container_type &dqdt ) const
+ {
+ for( size_t i=0 ; i<p.size() ; ++i ) dqdt[i] = p[i] / m_masses[i];
+ }
+};
+</pre>
+
+ 
+
+ 
+
+ 
+
+
+
+<pre class="programlisting">struct solar_system
+{
+ const mass_type &m_masses;
+
+ solar_system( const mass_type &masses ) : m_masses( masses ) { }
+
+ void operator()( const container_type &q , container_type &dpdt ) const
+ {
+ const size_t n = q.size();
+ fill( dpdt.begin() , dpdt.end() , point_type( 0.0 , 0.0 , 0.0 ) );
+ for( size_t i=0 ; i<n ; ++i )
+ {
+ for( size_t j=i+1 ; j<n ; ++j )
+ {
+ point_type diff = q[j] - q[i];
+ double d = abs( diff );
+ diff *= ( gravitational_constant / d / d / d * m_masses[i] * m_masses[j] ) ;
+ dpdt[i] += diff ;
+ dpdt[j] -= diff ;
+ }
+ }
+ }
+};
+</pre>
+
+ 
+
 

 In general a three body-system is chaotic, hence we can not expect that
 arbitray initial conditions of the system will lead to a dynamic which
 is comparable with the solar system. That is we have to define proper initial
- conditions.
- 
-
- show the code
+ conditions, which are taken from the book of Hairer, Wannier, Lubich.
 

 Now, we use the rk stepper to integrate the solar system. To visualize
@@ -495,10 +575,35 @@
 The output can be piped directly into gnuplot and a very nice visualization
 of the motion appears.
 
-</div>

- usage of the steppers
- 
+ 
+
+
+
+<pre class="programlisting">typedef hamiltonian_stepper_euler< container_type > stepper_type;
+state_type state = make_pair( q , p );
+for( size_t c = 0 ; c<2000 ; ++c )
+{
+ clog << t << tab << energy( state.first , state.second , masses ) << tab;
+ clog << center_of_mass( state.first , masses ) << tab;
+ clog << center_of_mass( state.second , masses ) << endl;
+
+ cout << t;
+ for( size_t i=0 ; i<n ; ++i ) cout << tab << state.first[i];
+ cout << endl;
+
+ for( size_t i=0 ; i<100 ; ++i,t+=dt )
+ stepper.do_step( make_pair( solar_system_momentum( masses ) ,
+ solar_system( masses ) ) ,
+ state , t , dt );
+ t += dt;
+}
+</pre>
+
+ 
+
+ 
+</div>
</div>
<div class="section" lang="en">
<div class="titlepage"><div><div><h3 class="title">

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-07-05 10:50:10 EDT (Mon, 05 Jul 2010)
@@ -70,9 +70,10 @@
 x(t) are calculated iteratively. The easiest algorithm is the Euler-Scheme,
 where starting at x(0) one finds x(dt) = x(0) + dt*f(x(0),0). Now one can use
 x(dt) and obtain x(2dt) in a similar way and so on. The Euler method is of
- order 1, that means the error at each step is ~ dt[superscript 2]. This is,
- of course, not very satisfying, which is why the Euler method is merely used
- for real life problems and serves just as illustrative example. In odeint, the following algorithms are implemented:
+ order 1, that means the error at each step is ~ dt^2. This is, of course, not
+ very satisfying, which is why the Euler method is merely used for real life
+ problems and serves just as illustrative example. In odeint,
+ the following algorithms are implemented:
 
<div class="table">
<a name="boost_numeric_odeint.overview.stepper_algorithms"></a>Table 1.1. Stepper Algorithms
@@ -245,7 +246,7 @@
</div>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
-<td align="left">Last revised: July 05, 2010 at 13:27:18 GMT</td>
+<td align="left">Last revised: July 05, 2010 at 14:50:52 GMT</td>
<td align="right"><div class="copyright-footer"></div></td>
</tr></table>
<hr>

Modified: sandbox/odeint/libs/numeric/odeint/doc/tutorial.qbk
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/doc/tutorial.qbk (original)
+++ sandbox/odeint/libs/numeric/odeint/doc/tutorial.qbk 2010-07-05 10:50:10 EDT (Mon, 05 Jul 2010)
@@ -157,23 +157,23 @@
where pi is the momenta of object i. The equations of motion can also be
derived from the Hamiltonian

-H = sum over i pi^2 / 2 + sum j V( qi , qj )
+H = sum over i pi^2 / 2 m_i + sum j V( qi , qj )

via dqi = dH / dpi, dpi = - dH / dq1. V(qi,qj) is the interaction
potential.

In time independent Hamiltonian system the energy is conserved and special
integration methods have to be applied in order to ensure energy
-conservation. The odeint library provides two stepper classes for Hamiltonian
-systems, which are separable and can be written in the form H = sum pi^2/2 +
-Hq.
+conservation. The odeint library provides classes for Hamiltonian
+systems, which are separable and can be written in the form H = sum pi^2/2 mi +
+Hq(q), where Hq(q) only depends on the coordinates.

hamiltonian_stepper_euler
hamiltonian_stepper_rk

Alltough this functional form might look a bit arbitrary it covers nearly all
classical mechanical systems with inertia and without dissipation, or where
-the equations of motion can be written in the form dqi=pi dpi=f(qi).
+the equations of motion can be written in the form dqi=mi pi dpi=f(qi).

[endsect]
@@ -185,35 +185,38 @@
as well as the velocity. Therefore, we use the operators from
<boost/operator.hpp>:

-show the code
+[import ../examples/point_type.hpp]
+[point_type]

-The next step is to define the state type and the system (derivative)
-function. As state type we use std::tr1::array and a state type represents all
-space coordinates q or all momenta coordinates p. As system function we have
-to provide f(q)
+The next step is to define a container type storing the values of q and p and
+to define systems (derivative) functions. As container type we use
+std::tr1::array and the state type is then simply

-show the code
+[import ../examples/solar_system.cpp]
+[state_type_definition]

-Note, that we have allready define the masses of all planets in the solar
-system.
+and represents all space coordinates q or all momenta coordinates p. As system
+function we have to provide f(p) = dq and f(q) = -dp, which acts only on p or q:
+
+
+[momentum_function]
+
+[coordinate_function]

In general a three body-system is chaotic, hence we can not expect that
arbitray initial conditions of the system will lead to a dynamic which is
comparable with the solar system. That is we have to define proper initial
-conditions.
-
-show the code
+conditions, which are taken from the book of Hairer, Wannier, Lubich.

Now, we use the rk stepper to integrate the solar system. To visualize the
motion we save the trajectory of each planet in a circular buffer. The output
can be piped directly into gnuplot and a very nice visualization of the motion
appears.

-[endsect]
-
+[integration_solar_system]

-usage of the steppers
+[endsect]

[endsect]

Modified: sandbox/odeint/libs/numeric/odeint/examples/point_type.hpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/point_type.hpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/point_type.hpp 2010-07-05 10:50:10 EDT (Mon, 05 Jul 2010)
@@ -21,6 +21,8 @@
//
// the point type
//
+
+//[ point_type
template< class T , size_t Dim >
class point :
 boost::additive1< point< T , Dim > ,
@@ -28,8 +30,11 @@
 boost::multiplicative2< point< T , Dim > , T
> > >
{
+ // ...
+//<-
public:

+
 const static size_t dim = Dim;
 typedef T value_type;
 typedef point< value_type , dim > point_type;
@@ -99,8 +104,13 @@
private:

 T m_val[dim];
+
+//->
};

+// define the scalar_prod(), norm() and abs()
+//]
+

//
// the - operator

Modified: sandbox/odeint/libs/numeric/odeint/examples/solar_system.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/solar_system.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/solar_system.cpp 2010-07-05 10:50:10 EDT (Mon, 05 Jul 2010)
@@ -29,25 +29,28 @@
// we simulate 5 planets and the sun
const size_t n = 6;

+//[ state_type_definition
typedef point< double , 3 > point_type;
-typedef std::tr1::array< point_type , n > state_type;
+typedef std::tr1::array< point_type , n > container_type;
typedef std::tr1::array< double , n > mass_type;
-typedef hamiltonian_stepper_euler< state_type > stepper_type;
+typedef pair< container_type , container_type > state_type;
+//]
+

const double gravitational_constant = 2.95912208286e-4;

-ostream& operator<<( ostream &out , const state_type &x )
+ostream& operator<<( ostream &out , const container_type &x )
{
- typedef state_type::value_type value_type;
+ typedef container_type::value_type value_type;
 copy( x.begin() , x.end() ,
 ostream_iterator< value_type >( out , "\n" ) );
 return out;
}

-point_type center_of_mass( const state_type &x , const mass_type &m )
+point_type center_of_mass( const container_type &x , const mass_type &m )
{
 double overall_mass = 0.0;
 point_type mean( 0.0 );
@@ -60,7 +63,7 @@
 return mean;
}

-point_type center( const state_type &x )
+point_type center( const container_type &x )
{
 point_type mean( 0.0 );
 for( size_t i=0 ; i<x.size() ; ++i ) mean += x[i];
@@ -69,7 +72,7 @@
}

-double energy( const state_type &q , const state_type &p , const mass_type &masses )
+double energy( const container_type &q , const container_type &p , const mass_type &masses )
{
 const size_t n = q.size();
 double en = 0.0;
@@ -85,13 +88,15 @@
 return en;
}

+
+//[ coordinate_function
struct solar_system
{
 const mass_type &m_masses;

 solar_system( const mass_type &masses ) : m_masses( masses ) { }

- void operator()( const state_type &q , state_type &dpdt ) const
+ void operator()( const container_type &q , container_type &dpdt ) const
 {
 const size_t n = q.size();
 fill( dpdt.begin() , dpdt.end() , point_type( 0.0 , 0.0 , 0.0 ) );
@@ -108,19 +113,21 @@
 }
 }
};
+//]

+//[ momentum_function
struct solar_system_momentum
{
 const mass_type &m_masses;

 solar_system_momentum( const mass_type &masses ) : m_masses( masses ) { }

- void operator()( const state_type &p , state_type &dqdt ) const
+ void operator()( const container_type &p , container_type &dqdt ) const
 {
 for( size_t i=0 ; i<p.size() ; ++i ) dqdt[i] = p[i] / m_masses[i];
 }
};
-
+//]

@@ -135,7 +142,7 @@
 1.0 / ( 1.3e8 ) // pluto
 }};

- state_type q = {{
+ container_type q = {{
 point_type( 0.0 , 0.0 , 0.0 ) , // sun
 point_type( -3.5023653 , -3.8169847 , -1.5507963 ) , // jupiter
 point_type( 9.0755314 , -3.0458353 , -1.6483708 ) , // saturn
@@ -144,7 +151,7 @@
 point_type( -15.5387357 , -25.2225594 , -3.1902382 ) // pluto
 }};

- state_type p = {{
+ container_type p = {{
 point_type( 0.0 , 0.0 , 0.0 ) , // sun
 point_type( 0.00565429 , -0.00412490 , -0.00190589 ) , // jupiter
 point_type( 0.00168318 , 0.00483525 , 0.00192462 ) , // saturn
@@ -162,8 +169,9 @@
 const double dt = 1.0;
 double t = 0.0;

-
- pair< state_type , state_type > state = make_pair( q , p );
+//[ integration_solar_system
+ typedef hamiltonian_stepper_euler< container_type > stepper_type;
+ state_type state = make_pair( q , p );
 for( size_t c = 0 ; c<2000 ; ++c )
 {
 clog << t << tab << energy( state.first , state.second , masses ) << tab;
@@ -180,6 +188,7 @@
 state , t , dt );
 t += dt;
 }
+//]

 return 0;
}

Next message: karsten.ahnert_at_[hidden]: "[Boost-commit] svn:boost r63653 - in sandbox/odeint/libs/numeric/odeint/doc: . html html/boost_numeric_odeint"
Previous message: anthony_at_[hidden]: "[Boost-commit] svn:boost r63651 - trunk/libs/thread/test"

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