Subject: [Boost-commit] svn:boost r57461 - in sandbox/odeint: boost/numeric boost/numeric/odeint libs/numeric/odeint/examples
From: karsten.ahnert_at_[hidden]
Date: 2009-11-07 14:05:47
Author: karsten
Date: 2009-11-07 14:05:44 EST (Sat, 07 Nov 2009)
New Revision: 57461
URL: http://svn.boost.org/trac/boost/changeset/57461
Log:
half_step stepper included and clean up examples
Added:
sandbox/odeint/libs/numeric/odeint/examples/lorenz_stepper.cpp (contents, props changed)
- copied, changed from r57443, /sandbox/odeint/libs/numeric/odeint/examples/lorenz.cpp
Removed:
sandbox/odeint/libs/numeric/odeint/examples/lorenz.cpp
sandbox/odeint/libs/numeric/odeint/examples/lorenz_array.cpp
sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp
Text files modified:
sandbox/odeint/boost/numeric/odeint.hpp | 3 +
sandbox/odeint/boost/numeric/odeint/euler.hpp | 45 +--------------------------
sandbox/odeint/boost/numeric/odeint/runge_kutta_4.hpp | 39 ++++++++++++++++--------
sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp | 19 ++++++-----
sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp | 8 +++-
sandbox/odeint/libs/numeric/odeint/examples/lorenz_stepper.cpp | 64 ++++++++++++++++++++++++++++-----------
6 files changed, 91 insertions(+), 87 deletions(-)
Modified: sandbox/odeint/boost/numeric/odeint.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint.hpp 2009-11-07 14:05:44 EST (Sat, 07 Nov 2009)
@@ -17,7 +17,10 @@
#include <boost/numeric/odeint/euler.hpp>
#include <boost/numeric/odeint/runge_kutta_4.hpp>
+#include <boost/numeric/odeint/stepper_half_step.hpp>
+
#include <boost/numeric/odeint/stepsize_controller_standard.hpp>
+
#include <boost/numeric/odeint/integrator.hpp>
#include <boost/numeric/odeint/integrator_constant_step.hpp>
Modified: sandbox/odeint/boost/numeric/odeint/euler.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/euler.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/euler.hpp 2009-11-07 14:05:44 EST (Sat, 07 Nov 2009)
@@ -44,6 +44,7 @@
typedef Container container_type;
typedef Resizer resizer_type;
typedef Time time_type;
+ typedef const unsigned short order_type;
typedef typename container_type::value_type value_type;
typedef typename container_type::iterator iterator;
@@ -72,7 +73,7 @@
// public interface
public:
- const unsigned int order() { return 1; }
+ order_type order() const { return 1; }
@@ -86,8 +87,6 @@
detail::it_algebra::increment( x.begin() , x.end() , dxdt.begin() , dt );
}
-
-
template< class DynamicalSystem >
void next_step( DynamicalSystem system ,
container_type &x ,
@@ -98,46 +97,6 @@
system( x , m_dxdt , t );
next_step( system , x , m_dxdt , t , dt );
}
-
-
-
- template< class DynamicalSystem >
- void next_step( DynamicalSystem system ,
- container_type &x ,
- const container_type &dxdt ,
- time_type t ,
- time_type dt ,
- container_type &xerr )
- {
- m_resizer.adjust_size( x , xerr );
-
- m_xtemp = x;
- time_type dt2 = 0.5 * dt;
-
- next_step( system , x , dxdt , t , dt );
- next_step( system , m_xtemp , dxdt , t , dt2 );
- next_step( system , m_xtemp , t+dt2 , dt2 );
-
- //xerr = x - m_xtemp
- detail::it_algebra::assign_diff(xerr.begin() ,
- xerr.end() ,
- x.begin() ,
- m_xtemp.begin() );
- }
-
-
-
- template< class DynamicalSystem >
- void next_step( DynamicalSystem system ,
- container_type &x ,
- time_type t ,
- time_type dt ,
- container_type &xerr )
- {
- m_resizer.check_size_and_resize( x , m_dxdt );
- system( x , m_dxdt , t );
- next_step( system , x , m_dxdt , t , dt , xerr );
- }
};
Modified: sandbox/odeint/boost/numeric/odeint/runge_kutta_4.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/runge_kutta_4.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/runge_kutta_4.hpp 2009-11-07 14:05:44 EST (Sat, 07 Nov 2009)
@@ -42,6 +42,7 @@
typedef Container container_type;
typedef Resizer resizer_type;
typedef Time time_type;
+ typedef const unsigned short order_type;
typedef typename container_type::value_type value_type;
typedef typename container_type::iterator iterator;
@@ -74,9 +75,12 @@
// public interface
public:
+ order_type order() const { return 4; }
+
template< class DynamicalSystem >
void next_step( DynamicalSystem system ,
container_type &x ,
+ const container_type &dxdt ,
time_type t ,
time_type dt )
{
@@ -84,30 +88,39 @@
const time_type val2 = time_type( 2.0 );
- m_resizer.adjust_size( x , m_dxdt );
- m_resizer.adjust_size( x , m_dxt );
- m_resizer.adjust_size( x , m_dxm );
- m_resizer.adjust_size( x , m_xt );
+ m_resizer.adjust_size( x , m_dxt );
+ m_resizer.adjust_size( x , m_dxm );
+ m_resizer.adjust_size( x , m_xt );
time_type dh = time_type( 0.5 ) * dt;
time_type th = t + dh;
- system( x , m_dxdt , t );
- //m_xt = x + dh*m_dxdt
- assign_sum( m_xt.begin() , m_xt.end() , x.begin() , m_dxdt.begin() , dh );
+ assign_sum( m_xt.begin() , m_xt.end() , x.begin() , dxdt.begin() , dh );
system( m_xt , m_dxt , th );
- //m_xt = x + dh*m_dxdt
assign_sum( m_xt.begin() , m_xt.end() , x.begin() , m_dxt.begin() , dh );
system( m_xt , m_dxm , th );
- //m_xt = x + dt*m_dxm ; m_dxm += m_dxt
- assign_sum_increment( m_xt.begin() , m_xt.end() , x.begin() , m_dxm.begin() , m_dxt.begin() , dt );
+ assign_sum_increment( m_xt.begin() , m_xt.end() , x.begin() ,
+ m_dxm.begin() , m_dxt.begin() , dt );
system( m_xt , m_dxt , value_type( t + dt ) );
- //x = dt/6 * ( m_dxdt + m_dxt + val2*m_dxm )
- increment_sum_sum( x.begin() , x.end() , m_dxdt.begin() , m_dxt.begin() , m_dxm.begin() , dt / time_type( 6.0 ) , val2
- );
+ increment_sum_sum( x.begin() , x.end() , dxdt.begin() ,
+ m_dxt.begin() , m_dxm.begin() ,
+ dt / time_type( 6.0 ) , val2 );
+ }
+
+
+
+ template< class DynamicalSystem >
+ void next_step( DynamicalSystem system ,
+ container_type &x ,
+ time_type t ,
+ time_type dt )
+ {
+ m_resizer.adjust_size( x , m_dxdt );
+ system( x , m_dxdt , t );
+ next_step( system , x , m_dxdt , t , dt );
}
Modified: sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp
==============================================================================
--- sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp (original)
+++ sandbox/odeint/boost/numeric/odeint/stepper_half_step.hpp 2009-11-07 14:05:44 EST (Sat, 07 Nov 2009)
@@ -13,8 +13,8 @@
copy at http://www.boost.org/LICENSE_1_0.txt)
*/
-#ifndef BOOST_NUMERIC_ODEINT_EULER_HPP
-#define BOOST_NUMERIC_ODEINT_EULER_HPP
+#ifndef BOOST_NUMERIC_ODEINT_STEPPER_HALF_STEP_HPP
+#define BOOST_NUMERIC_ODEINT_STEPPER_HALF_STEP_HPP
#include <boost/concept_check.hpp>
@@ -41,6 +41,7 @@
typedef typename Stepper::container_type container_type;
typedef typename Stepper::resizer_type resizer_type;
typedef typename Stepper::time_type time_type;
+ typedef typename Stepper::order_type order_type;
typedef typename container_type::value_type value_type;
typedef typename container_type::iterator iterator;
@@ -59,7 +60,10 @@
stepper_type m_stepper;
- const unsigned int order() { return m_stepper.order(); }
+ // public interface
+ public:
+
+ order_type order() const { return m_stepper.order(); }
@@ -84,8 +88,6 @@
m_stepper.next_step( system , x , t , dt );
}
- /*
-
template< class DynamicalSystem >
void next_step( DynamicalSystem system ,
container_type &x ,
@@ -119,12 +121,11 @@
time_type dt ,
container_type &xerr )
{
- m_resizer.check_size_and_resize( x , m_dxdt );
+ m_resizer.adjust_size( x , m_dxdt );
system( x , m_dxdt , t );
next_step( system , x , m_dxdt , t , dt , xerr );
}
- */
- }
+ };
@@ -133,4 +134,4 @@
} // namespace boost
-#endif // BOOST_NUMERIC_ODEINT_EULER_HPP
+#endif // BOOST_NUMERIC_ODEINT_STEPPER_HALF_STEP_HPP
Deleted: sandbox/odeint/libs/numeric/odeint/examples/lorenz.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz.cpp 2009-11-07 14:05:44 EST (Sat, 07 Nov 2009)
+++ (empty file)
@@ -1,76 +0,0 @@
-/* Boost numeric/odeint/examples/lorenz.cpp
-
- Copyright 2009 Karsten Ahnert
-
- Shows the usage of odeint, and integrates the Lorenz equations,
- a deterministic chaotic system
-
- dx/dt = sigma * ( x - y)
- dy/dt = R*x - y - x*z
- dz/dt = x*y - b*z
-
- mit sigma = 10, r=28, b = 8/3
-
- 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)
-*/
-
-#include <iostream>
-#include <vector>
-#include <list>
-#include <tr1/array>
-
-#include <boost/numeric/odeint.hpp>
-
-#define tab "\t"
-
-using namespace std;
-using namespace std::tr1;
-using namespace boost::numeric::odeint;
-
-
-
-
-const double sigma = 10.0;
-const double R = 28.0;
-const double b = 8.0 / 3.0;
-
-const double dt = 0.01;
-const size_t olen = 10000;
-
-typedef vector<double> state_type;
-
-void lorenz( state_type &x , state_type &dxdt , double t )
-{
- dxdt[0] = sigma * ( x[1] - x[0] );
- dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
- dxdt[2] = x[0]*x[1] - b * x[2];
-}
-
-int main( int argc , char **argv )
-{
-
- state_type x(3);
- x[0] = 1.0;
- x[1] = 0.0;
- x[2] = 0.0;
-
- ode_step_euler< state_type > euler;
-
- double t = 0.0;
- for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
- {
- cout << t << tab << x[0] << tab << x[1] << tab << x[2] << endl;
- euler.next_step( lorenz , x , t , dt );
- }
-
- return 0;
-}
-
-
-
-/*
- Compile with
- g++ -Wall -I$BOOST_ROOT -I../../../../ lorenz.cpp
-*/
Deleted: sandbox/odeint/libs/numeric/odeint/examples/lorenz_array.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_array.cpp 2009-11-07 14:05:44 EST (Sat, 07 Nov 2009)
+++ (empty file)
@@ -1,76 +0,0 @@
-/* Boost numeric/odeint/examples/lorenz_array.cpp
-
- Copyright 2009 Karsten Ahnert
-
- Shows, the usage of odeint, and integrates the Lorenz equations,
- a deterministic chaotic system
-
- dx/dt = sigma * ( x - y)
- dy/dt = R*x - y - x*z
- dz/dt = x*y - b*z
-
- mit sigma = 10, r=28, b = 8/3
-
- 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)
-*/
-
-#include <iostream>
-#include <vector>
-#include <list>
-#include <tr1/array>
-
-#include <boost/numeric/odeint.hpp>
-
-#define tab "\t"
-
-using namespace std;
-using namespace std::tr1;
-using namespace boost::numeric::odeint;
-
-const double sigma = 10.0;
-const double R = 28.0;
-const double b = 8.0 / 3.0;
-
-const double dt = 0.01;
-const size_t olen = 10000;
-
-typedef array<double, 3> state_type;
-
-void lorenz( state_type &x , state_type &dxdt , double t )
-{
- dxdt[0] = sigma * ( x[1] - x[0] );
- dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
- dxdt[2] = x[0]*x[1] - b * x[2];
-}
-
-int main( int argc , char **argv )
-{
- state_type x , x2;
- x[0] = 1.0;
- x[1] = 0.0;
- x[2] = 20.0;
- x2 = x;
-
- ode_step_euler< state_type > euler;
- ode_step_runge_kutta_4< state_type > rk4;
-
- double t = 0.0;
- for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
- {
- cout << t << tab << x[0] << tab << x[1] << tab << x[2] << tab;
- cout << x2[0] << tab << x2[1] << tab << x2[2] << endl;
- euler.next_step( lorenz , x , t , dt );
- rk4.next_step( lorenz , x2 , t , dt );
- }
-
- return 0;
-}
-
-
-
-/*
- Compile with
- g++ -Wall -I$BOOST_ROOT -I../../../../ lorenz_array.cpp
-*/
Modified: sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_controlled.cpp 2009-11-07 14:05:44 EST (Sat, 07 Nov 2009)
@@ -1,6 +1,7 @@
/* Boost numeric/odeint/examples/lorenz_controlled.cpp
Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
Shows the usage of odeint by integrating the Lorenz equations,
a deterministic chaotic system
@@ -59,9 +60,10 @@
x[1] = 0.0;
x[2] = 20.0;
- ode_step_euler< state_type > euler;
+ ode_step_half_step< ode_step_euler< state_type > > euler;
+// ode_step_euler< state_type > euler;
step_controller_standard< state_type, double >
- controller(eps_abs, eps_rel, 1.0, 1.0);
+ controller( eps_abs , eps_rel, 1.0, 1.0);
cout.precision(5);
cout.setf(ios::fixed,ios::floatfield);
@@ -75,5 +77,5 @@
/*
Compile with
- g++ -Wall -I$BOOST_ROOT -I../../../../ lorenz_array.cpp
+ g++ -Wall -O3 -I$BOOST_ROOT -I../../../../ lorenz_controlled.cpp
*/
Deleted: sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp
==============================================================================
--- sandbox/odeint/libs/numeric/odeint/examples/lorenz_integrator.cpp 2009-11-07 14:05:44 EST (Sat, 07 Nov 2009)
+++ (empty file)
@@ -1,85 +0,0 @@
-/* Boost numeric/odeint/examples/lorenz_integrator.cpp
-
- Copyright 2009 Karsten Ahnert
- Copyright 2009 Mario Mulansky
-
- Shows the usage of odeint integrator by integrating the Lorenz equations,
- a deterministic chaotic system
-
- dx/dt = sigma * ( x - y)
- dy/dt = R*x - y - x*z
- dz/dt = x*y - b*z
-
- mit sigma = 10, r=28, b = 8/3
-
- 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)
-*/
-
-#include <iostream>
-#include <vector>
-#include <list>
-#include <tr1/array>
-
-#include <boost/numeric/odeint.hpp>
-
-#define tab "\t"
-
-using namespace std;
-using namespace std::tr1;
-using namespace boost::numeric::odeint;
-
-const double sigma = 10.0;
-const double R = 28.0;
-const double b = 8.0 / 3.0;
-
-const double eps_abs = 1E-3;
-const double eps_rel = 1E-3;
-
-const size_t time_points = 201;
-
-typedef array<double, 3> state_type;
-
-void lorenz( state_type &x , state_type &dxdt , double t )
-{
- dxdt[0] = sigma * ( x[1] - x[0] );
- dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
- dxdt[2] = x[0]*x[1] - b * x[2];
-}
-
-int main( int argc , char **argv )
-{
- state_type x;
- x[0] = 1.0;
- x[1] = 0.0;
- x[2] = 20.0;
-
- vector<state_type> x_t_vec(time_points);
- vector<double> times(time_points);
- for( size_t i=0; i<time_points; i++ ) {
- times[i] = 0.1*i;
- }
-
- ode_step_euler< state_type > euler;
- size_t steps = integrate( euler, lorenz, x, times, x_t_vec);
-
- clog << "Steps: " << steps << endl;
-
- cout.precision(5);
- cout.setf(ios::fixed,ios::floatfield);
-
-
- for( size_t i=0; i<time_points; 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;
- }
-
- return 0;
-}
-
-/*
- Compile with
- g++ -Wall -I$BOOST_ROOT -I../../../../ lorenz_integrator.cpp
-*/
Copied: sandbox/odeint/libs/numeric/odeint/examples/lorenz_stepper.cpp (from r57443, /sandbox/odeint/libs/numeric/odeint/examples/lorenz.cpp)
==============================================================================
--- /sandbox/odeint/libs/numeric/odeint/examples/lorenz.cpp (original)
+++ sandbox/odeint/libs/numeric/odeint/examples/lorenz_stepper.cpp 2009-11-07 14:05:44 EST (Sat, 07 Nov 2009)
@@ -1,6 +1,7 @@
-/* Boost numeric/odeint/examples/lorenz.cpp
+/* Boost numeric/odeint/examples/lorenz_stepper.cpp
Copyright 2009 Karsten Ahnert
+ Copyright 2009 Mario Mulansky
Shows the usage of odeint, and integrates the Lorenz equations,
a deterministic chaotic system
@@ -9,7 +10,10 @@
dy/dt = R*x - y - x*z
dz/dt = x*y - b*z
- mit sigma = 10, r=28, b = 8/3
+ with sigma = 10, r=28, b = 8/3
+
+ Furthermore, the usage of std::tr1::array and std::vector in odeint is
+ shown and the performance of both containers is compared.
Distributed under the Boost Software License, Version 1.0.
(See accompanying file LICENSE_1_0.txt or
@@ -26,44 +30,66 @@
#define tab "\t"
using namespace std;
-using namespace std::tr1;
using namespace boost::numeric::odeint;
+typedef std::vector< double > state_type1;
+typedef std::tr1::array< double , 3 > state_type2;
const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;
-const double dt = 0.01;
-const size_t olen = 10000;
-
-typedef vector<double> state_type;
+void lorenz1( state_type1 &x , state_type1 &dxdt , double t )
+{
+ dxdt[0] = sigma * ( x[1] - x[0] );
+ dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+ dxdt[2] = x[0]*x[1] - b * x[2];
+}
-void lorenz( state_type &x , state_type &dxdt , double t )
+void lorenz2( state_type2 &x , state_type2 &dxdt , double t )
{
dxdt[0] = sigma * ( x[1] - x[0] );
dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
dxdt[2] = x[0]*x[1] - b * x[2];
}
+
+
+
int main( int argc , char **argv )
{
+ const double dt = 0.01;
+ const size_t olen = 1000000000;
- state_type x(3);
- x[0] = 1.0;
- x[1] = 0.0;
- x[2] = 0.0;
+ state_type1 x1(3);
+ x1[0] = 1.0;
+ x1[1] = 0.0;
+ x1[2] = 0.0;
+ state_type2 x2 = {{ 1.0 , 0.0 , 0.0 }};
- ode_step_euler< state_type > euler;
+ ode_step_runge_kutta_4< state_type1 > stepper1;
+ ode_step_runge_kutta_4< state_type2 > stepper2;
- double t = 0.0;
+ clock_t start , end;
+ double t;
+
+ start= clock();
+ t = 0.0;
for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
- {
- cout << t << tab << x[0] << tab << x[1] << tab << x[2] << endl;
- euler.next_step( lorenz , x , t , dt );
- }
+ stepper1.next_step( lorenz1 , x1 , t , dt );
+ end = clock();
+ cout << "vector : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
+
+
+ start= clock();
+ t = 0.0;
+ for( size_t oi=0 ; oi<olen ; ++oi,t+=dt )
+ stepper2.next_step( lorenz2 , x2 , t , dt );
+ end = clock();
+ cout << "array : " << double ( end - start ) / double( CLOCKS_PER_SEC ) << endl;
+
return 0;
}
@@ -72,5 +98,5 @@
/*
Compile with
- g++ -Wall -I$BOOST_ROOT -I../../../../ lorenz.cpp
+ g++ -Wall -O3 -I$BOOST_ROOT -I../../../../ lorenz_stepper.cpp
*/